PeriDyno 1.0.0
Loading...
Searching...
No Matches
ImplicitISPH.cu
Go to the documentation of this file.
1#include "ImplicitISPH.h"
2
3#include "SummationDensity.h"
4
5namespace dyno
6{
7 template<typename TDataType>
8 ImplicitISPH<TDataType>::ImplicitISPH()
9 : ParticleApproximation<TDataType>()
10 {
11 //this->varIterationNumber()->setValue(3);
12 this->varKappa()->setValue(200);
13 this->varRestDensity()->setValue(Real(1000));
14
15 mSummation = std::make_shared<SummationDensity<TDataType>>();
16
17 this->inSmoothingLength()->connect(mSummation->inSmoothingLength());
18 this->inSamplingDistance()->connect(mSummation->inSamplingDistance());
19 this->inPosition()->connect(mSummation->inPosition());
20 this->inNeighborIds()->connect(mSummation->inNeighborIds());
21
22 mSummation->outDensity()->connect(this->outDensity());
23 }
24
25 template<typename TDataType>
26 ImplicitISPH<TDataType>::~ImplicitISPH()
27 {
28 mSourceTerm.clear();
29 mDii.clear();
30 mAii.clear();
31 mAnPn.clear();
32 mSumDijPj.clear();
33 mPressrue.clear();
34 mOldPressrue.clear();
35 m_Residual.clear();
36 mPredictDensity.clear();
37 mDensityAdv.clear();
38
39 }
40
41
42 /*
43 *@Equation: Source_i = rho_0 - (rho_i + dt * div(vi))
44 */
45 template <typename Real, typename Coord, typename Kernel>
46 __global__ void IISPH_SourceTermCompute(
47 DArray<Real> Sources,
48 DArray<Coord> posArr,
49 DArray<Coord> veloArr,
50 DArray<Real> rhoArr,
51 DArrayList<int> neighbors,
52 Real dt,
53 Real rest_density,
54 Real mass,
55 Real smoothingLength,
56 Kernel gradient,
57 Real scale
58 ) {
59
60 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
61 if (pId >= Sources.size()) return;
62
63 Coord pos_i = posArr[pId];
64
65 Real div_i = 0.0f;
66
67 List<int>& list_i = neighbors[pId];
68 int nbSize = list_i.size();
69
70 for (int ne = 0; ne < nbSize; ne++)
71 {
72 int j = list_i[ne];
73 Real r = (pos_i - posArr[j]).norm();
74
75 if (r > EPSILON)
76 {
77 Coord g = gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
78 div_i += mass * g.dot(veloArr[pId] - veloArr[j]);
79 }
80 }
81
82 Sources[pId] = rest_density - (rhoArr[pId] + dt * div_i);
83
84 if (Sources[pId] > 0.0f) Sources[pId] = 0.0f;
85
86 }
87
88 /*
89 *@Equation: Coord D_ii = - dt^2 * sum_j { mass/(rho_i)^2 * Grad_Wij }
90 */
91 template <typename Real, typename Coord, typename Kernel>
92 __global__ void IISPH_DiiCoefficientOfPpeCompute(
93 DArray<Coord> D,
94 DArray<Coord> posArr,
95 DArray<Real> rhoArr,
96 DArrayList<int> neighbors,
97 Real dt,
98 Real rest_density,
99 Real mass,
100 Real smoothingLength,
101 Kernel gradient,
102 Real scale
103 ) {
104
105 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
106 if (pId >= D.size()) return;
107
108 Coord pos_i = posArr[pId];
109
110 Coord sum(0.0f);
111
112 List<int>& list_i = neighbors[pId];
113 int nbSize = list_i.size();
114
115 for (int ne = 0; ne < nbSize; ne++)
116 {
117 int j = list_i[ne];
118 Real r = (pos_i - posArr[j]).norm();
119
120 if (r > EPSILON)
121 {
122 Coord g = gradient(r, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r);
123 sum -= mass / (rhoArr[pId] * rhoArr[pId]) * g;
124 }
125 }
126
127 D[pId] = sum * dt * dt;
128 }
129
130
131
132 /*
133 *@Equation: Coord A_ii = - sum_j { mass (D_ii - D_ji) Grad_Wij };
134 * D_ij = - dt^2 * mass / rho_i^2 * Grad_Wij;
135 */
136 template <typename Real, typename Coord, typename Kernel>
137 __global__ void IISPH_AiiCoefficientOfPpeCompute(
138 DArray<Real> Aii,
139 DArray<Coord> Dii,
140 DArray<Coord> posArr,
141 DArray<Real> rhoArr,
142 DArrayList<int> neighbors,
143 Real dt,
144 Real rest_density,
145 Real mass,
146 Real smoothingLength,
147 Kernel gradient,
148 Real scale
149 ) {
150 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
151 if (pId >= Aii.size()) return;
152
153 Coord pos_i = posArr[pId];
154
155 Real sum = (0.0f);
156
157 List<int>& list_i = neighbors[pId];
158 int nbSize = list_i.size();
159
160 Coord dji(0.0f);
161
162 for (int ne = 0; ne < nbSize; ne++)
163 {
164 int j = list_i[ne];
165 Real r = (pos_i - posArr[j]).norm();
166
167 if (r > EPSILON)
168 {
169 Coord g_ij = gradient(r, smoothingLength, scale) * (pos_i - posArr[j] ) * (1.0f / r);
170 Coord g_ji = gradient(r, smoothingLength, scale) * (posArr[j] - pos_i) * (1.0f / r);
171 dji = -1.0 * dt * dt * mass * g_ji / (rhoArr[pId] * rhoArr[pId]); //grad_W_ij ???(-1)
172 sum += mass * g_ij.dot(Dii[pId] - dji); //grad_W_ji ???
173 }
174 }
175
176 Aii[pId] = sum;
177 }
178
179
180 /*
181 *@Equation: Coord DijPj = D_ij * Pj;
182 * D_ij = - dt^2 * mass * Grad_Wij * Pj / rho_i^2 ;
183 */
184 template <typename Real, typename Coord, typename Kernel>
185 __global__ void IISPH_Sum_DijDotPjInPPE(
186 DArray<Coord> SumDijPj,
187 DArray<Real> pressures,
188 DArray<Coord> posArr,
189 DArray<Real> rhoArr,
190 DArrayList<int> neighbors,
191 Real dt,
192 Real rest_density,
193 Real mass,
194 Real smoothingLength,
195 Kernel gradient,
196 Real scale
197 ) {
198 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
199 if (pId >= SumDijPj.size()) return;
200
201 Coord pos_i = posArr[pId];
202
203 List<int>& list_i = neighbors[pId];
204 int nbSize = list_i.size();
205
206 Coord sum(0.0f);
207
208 for (int ne = 0; ne < nbSize; ne++)
209 {
210 int j = list_i[ne];
211 Real r_ij = (pos_i - posArr[j]).norm();
212 if (r_ij > EPSILON)
213 {
214 Coord g_ij = gradient(r_ij, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r_ij);
215 sum -= dt * dt * mass * g_ij * pressures[j] / (rhoArr[j] * rhoArr[j]);
216 }
217 }
218 SumDijPj[pId] = sum;
219 }
220
221 /*
222 *@Equation: Coord AnPn = sum_n {An * Pn} - Aij;
223 * AnPn = Sum_j { mass * {(Sum_j Dij * Pj) - D_jj * Pj - Sum_(k!=i) D_jk * Pj } Grad_Wij};
224 */
225 template <typename Real, typename Coord, typename Kernel>
226 __global__ void IISPH_AnDotPnInPPE(
227 DArray<Real> AnPn,
228 DArray<Real> pressures,
229 DArray<Coord> Dii,
230 DArray<Coord> SumDijPj,
231 DArray<Coord> posArr,
232 DArray<Real> rhoArr,
233 DArrayList<int> neighbors,
234 Real dt,
235 Real rest_density,
236 Real mass,
237 Real smoothingLength,
238 Kernel gradient,
239 Real scale
240 ) {
241
242 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
243 if (pId >= AnPn.size()) return;
244
245 Coord pos_i = posArr[pId];
246
247 List<int>& list_i = neighbors[pId];
248 int nbSize = list_i.size();
249
250 Coord d_ij(0.0f);
251 Coord d_jk(0.0f);
252
253 Real sum = 0.0f;
254
255 Coord sum_DjkPk(0.0f);
256
257 Coord d_ji(0.0f);
258
259 for (int ne = 0; ne < nbSize; ne++)
260 {
261 int j = list_i[ne];
262
263 Real r_ij = (pos_i - posArr[j]).norm();
264
265 if (r_ij > EPSILON)
266 {
267 Coord g_ij = gradient(r_ij, smoothingLength, scale) * (pos_i - posArr[j]) * (1.0f / r_ij);
268
269 //Coord g_ji = gradient(r_ij, smoothingLength, scale) * (posArr[j] - pos_i) * (1.0f / r_ij);
270
271 //Coord d_ji = -dt * dt * mass * g_ji / (rhoArr[pId] * rhoArr[pId]);
272
273 //sum += mass * g_ij.dot(SumDijPj[pId] - Dii[j] * pressures[j] - (SumDijPj[j] - d_ji * pressures[pId] ));
274
275 //Coord sumDjkPk = SumDijPj[j];
276 Coord sumDjkPk(0);
277
278 List<int>& list_j = neighbors[j];
279 for (int nj = 0; nj < list_j.size(); nj++)
280 {
281 int k = list_i[nj];
282 Real r_jk = (posArr[j] - posArr[k]).norm();
283 if (r_jk > EPSILON && k != pId)
284 {
285 Coord g_jk = gradient(r_jk, smoothingLength, scale) * (posArr[j] - posArr[k]) * (1.0f / r_jk);
286 d_jk = (-dt * dt * mass / (rhoArr[k] * rhoArr[k])) * g_jk;
287
288 sumDjkPk += d_jk * pressures[k];
289 }
290 }
291
292 sum += mass * g_ij.dot(SumDijPj[pId] - Dii[j] * pressures[j] - sumDjkPk);
293 }
294
295 }
296 AnPn[pId] = sum;
297 }
298
299 /*
300 *@Brief: Relaxed Jacobi Method.
301 *@Equation: pressure_new[i] = Omega * pressure_old[i] + (1 - Omega) * (source[pId] - AnPn) / Aii
302 *
303 */
304 template <typename Real>
305 __global__ void IISPH_PressureJacobiCompute(
306 DArray<Real> pressures,
307 DArray<Real> oldPressures,
308 DArray<Real> AnPn,
309 DArray<Real> Source,
310 DArray<Real> Aii,
311 DArray<Real> Residual,
312 Real omega
313 ) {
314
315 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
316 if (pId >= pressures.size()) return;
317
318
319 if (abs(Aii[pId]) > EPSILON)
320 {
321 pressures[pId] = oldPressures[pId] * (1.0 - omega) + omega * (Source[pId] - AnPn[pId]) / Aii[pId];
322 }
323 else {
324 pressures[pId] = 0.0f;
325 }
326
327 /*
328 *@note Clamping negative pressures will lead to incorrect estimation of the residual in the PPE.
329 */
330 //if (pressures[pId] < 0.0f)
331 // pressures[pId] = 0.0f;
332
333 Residual[pId] = abs(Source[pId] - AnPn[pId] - Aii[pId]* pressures[pId]);
334
335
336 }
337
338
339 template<typename TDataType>
340 void ImplicitISPH<TDataType>::compute()
341 {
342 int num = this->inPosition()->size();
343
344 if (this->outDensity()->size() != this->inPosition()->size())
345 this->outDensity()->resize(this->inPosition()->size());
346
347 if (mPredictDensity.size() != this->inPosition()->size())
348 mPredictDensity.resize(this->inPosition()->size());
349
350 if (mDensityAdv.size() != this->inPosition()->size())
351 mDensityAdv.resize(this->inPosition()->size());
352
353 if (mPressrue.size() != this->inPosition()->size())
354 mPressrue.resize(this->inPosition()->size());
355
356 if (mSourceTerm.size() != this->inPosition()->size())
357 mSourceTerm.resize(this->inPosition()->size());
358
359 if (mDii.size() != this->inPosition()->size())
360 mDii.resize(this->inPosition()->size());
361
362 if (mAii.size() != this->inPosition()->size())
363 mAii.resize(this->inPosition()->size());
364
365 if (mAnPn.size() != this->inPosition()->size())
366 mAnPn.resize(this->inPosition()->size());
367
368 if (mPressrue.size() != this->inPosition()->size())
369 mPressrue.resize(this->inPosition()->size());
370
371 if (mSumDijPj.size() != this->inPosition()->size())
372 mSumDijPj.resize(this->inPosition()->size());
373
374 if (m_Residual.size() != this->inPosition()->size())
375 m_Residual.resize(this->inPosition()->size());
376
377 if (mOldPressrue.size() != this->inPosition()->size())
378 mOldPressrue.resize(this->inPosition()->size());
379
380
381 mSummation->varRestDensity()->setValue(this->varRestDensity()->getValue());
382 mSummation->varKernelType()->setCurrentKey(this->varKernelType()->currentKey());
383 mSummation->update();
384
385 this->PreIterationCompute();
386
387 int it = 0;
388
389 int itNum = this->varIterationNumber()->getValue();
390
391 mOldPressrue.reset();
392 mPressrue.reset();
393
394// Real error_max = takeOneIteration();
395// Real error = error_max;
396
397 while ( (it++ < itNum))
398 {
399 takeOneIteration();
400 //std::cout << "Residual: " << error / error_max * 100.0f << "%" << std::endl;
401 }
402
403 updateVelocity();
404 }
405
406
407 template<typename TDataType>
408 void ImplicitISPH<TDataType>::PreIterationCompute()
409 {
410 Real dt = this->inTimeStep()->getData();
411 int num = this->inPosition()->size();
412 Real rho_0 = this->varRestDensity()->getValue();
413
414 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
415 IISPH_DiiCoefficientOfPpeCompute,
416 mDii,
417 this->inPosition()->getData(),
418 mSummation->outDensity()->getData(),
419 this->inNeighborIds()->getData(),
420 this->inTimeStep()->getValue(),
421 rho_0,
422 mSummation->getParticleMass(),
423 this->inSmoothingLength()->getValue()
424 );
425
426 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
427 IISPH_SourceTermCompute,
428 mSourceTerm,
429 this->inPosition()->getData(),
430 this->inVelocity()->getData(),
431 mSummation->outDensity()->getData(),
432 this->inNeighborIds()->getData(),
433 this->inTimeStep()->getValue(),
434 rho_0,
435 mSummation->getParticleMass(),
436 this->inSmoothingLength()->getValue()
437 );
438
439 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
440 IISPH_AiiCoefficientOfPpeCompute,
441 mAii,
442 mDii,
443 this->inPosition()->getData(),
444 mSummation->outDensity()->getData(),
445 this->inNeighborIds()->getData(),
446 this->inTimeStep()->getValue(),
447 rho_0,
448 mSummation->getParticleMass(),
449 this->inSmoothingLength()->getValue()
450 );
451
452
453 }
454
455
456 template<typename TDataType>
457 typename TDataType::Real ImplicitISPH<TDataType>::takeOneIteration()
458 {
459 Real dt = this->inTimeStep()->getData();
460 int num = this->inPosition()->size();
461 Real rho_0 = this->varRestDensity()->getValue();
462
463 mOldPressrue.assign(mPressrue);
464
465// cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
466// IISPH_Sum_DijDotPjInPPE,
467// mSumDijPj,
468// mPressrue,
469// this->inPosition()->getData(),
470// mSummation->outDensity()->getData(),
471// this->inNeighborIds()->getData(),
472// this->inTimeStep()->getValue(),
473// rho_0,
474// mSummation->getParticleMass(),
475// this->inSmoothingLength()->getValue()
476// );
477
478 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
479 IISPH_AnDotPnInPPE,
480 mAnPn,
481 mPressrue,
482 mDii,
483 mSumDijPj,
484 this->inPosition()->getData(),
485 mSummation->outDensity()->getData(),
486 this->inNeighborIds()->getData(),
487 this->inTimeStep()->getValue(),
488 rho_0,
489 mSummation->getParticleMass(),
490 this->inSmoothingLength()->getValue()
491 );
492
493
494 cuExecute(num, IISPH_PressureJacobiCompute,
495 mPressrue,
496 mOldPressrue,
497 mAnPn,
498 mSourceTerm,
499 mAii,
500 m_Residual,
501 this->varRelaxedOmega()->getValue()
502 );
503
504// auto m_reduce = Reduction<Real>::Create(num);
505// Real error = m_reduce->average(m_Residual.begin(), num) / num;
506 return 0;
507
508 }
509
510
511 template <typename Real, typename Coord, typename Kernel>
512 __global__ void IISPH_UpdateVelocity(
513 DArray<Coord> velocities,
514 DArray<Real> pressures,
515 DArray<Coord> posArr,
516 DArray<Real> rhoArr,
517 DArrayList<int> neighbors,
518 Real dt,
519 Real rest_density,
520 Real mass,
521 Real smoothingLength,
522 Kernel gradient,
523 Real scale
524 ) {
525 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
526 if (pId >= velocities.size()) return;
527
528 Coord Force_i(0.0f);
529
530 List<int>& list_i = neighbors[pId];
531 int nbSize = list_i.size();
532
533 for (int ne = 0; ne < nbSize; ne++)
534 {
535 int j = list_i[ne];
536 Real r_ij = (posArr[pId] - posArr[j]).norm();
537 if (r_ij > EPSILON)
538 {
539 Coord g_ij = gradient(r_ij, smoothingLength, scale) * (posArr[pId] - posArr[j]) * (1.0f / r_ij);
540 Force_i -= mass * mass * (pressures[pId] / (rhoArr[pId] * rhoArr[pId]) + pressures[j] / (rhoArr[j] * rhoArr[j])) * g_ij;
541 }
542 }
543
544 velocities[pId] += dt * Force_i / mass;
545 }
546
547 template<typename TDataType>
548 void ImplicitISPH<TDataType>::updateVelocity()
549 {
550 int num = this->inPosition()->size();
551 Real dt = this->inTimeStep()->getData();
552 Real rho_0 = this->varRestDensity()->getValue();
553
554 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
555 IISPH_UpdateVelocity,
556 this->inVelocity()->getData(),
557 mPressrue,
558 this->inPosition()->getData(),
559 mSummation->outDensity()->getData(),
560 this->inNeighborIds()->getData(),
561 this->inTimeStep()->getValue(),
562 rho_0,
563 mSummation->getParticleMass(),
564 this->inSmoothingLength()->getValue()
565 );
566
567
568 }
569
570 DEFINE_CLASS(ImplicitISPH);
571}
Generated by doxygen 1.13.2