PeriDyno 1.0.0
Loading...
Searching...
No Matches
SimpleVelocityConstraint.cu
Go to the documentation of this file.
1#include <cuda_runtime.h>
2#include "SimpleVelocityConstraint.h"
3#include <string>
4#include "Algorithm/Function2Pt.h"
5
6
7namespace dyno
8{
9 //F-norms
10 template<typename Real>
11 __device__ Real VB_FNorm(const Real a, const Real b, const Real c)
12 {
13 Real p = a + b + c;
14 return sqrt(p * p / 2);
15 }
16
17 //CrossModel Viscosity Coefficient
18 template<typename Real>
19 __device__ Real VB_Viscosity(const Real Viscosity_h, const Real Viscosity_l, const Real StrainRate, const Real CrossModel_K, const Real CrossModel_n)
20 {
21 Real p = CrossModel_K * StrainRate;
22 p = pow(p, CrossModel_n);
23 return Viscosity_h + (Viscosity_l - Viscosity_h) / (1 + p);
24 }
25
26 __device__ inline float kernWeight(const float r, const float h)
27 {
28 const float q = r / h;
29 if (q > 1.0f) return 0.0f;
30 else {
31 const float d = 1.0f - q;
32 const float hh = h * h;
33 return (1.0 - pow(q, 4.0f));
34 }
35 }
36
37
38 __device__ inline float kernWR(const float r, const float h)
39 {
40 float w = kernWeight(r, h);
41 const float q = r / h;
42 if (q < 0.4f)
43 {
44 return w / (0.4f * h);
45 }
46 return w / r;
47 }
48
49 __device__ inline float kernWRR(const float r, const float h)
50 {
51 float w = kernWeight(r, h);
52 const float q = r / h;
53 if (q < 0.4f)
54 {
55 return w / (0.16f * h * h);
56 }
57 return w / r / r;
58 }
59
60
61 template <typename Real, typename Coord>
62 __global__ void SIMPLE_ComputeAlpha
63 (
64 DArray<Real> alpha,
65 DArray<Coord> position,
66 DArray<Attribute> attribute,
67 DArrayList<int> neighbors,
68 Real smoothingLength
69 )
70 {
71 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
72 if (pId >= position.size()) return;
73 if (!attribute[pId].isFluid()) return;
74
75 Coord pos_i = position[pId];
76 Real alpha_i = 0.0f;
77
78 List<int>& list_i = neighbors[pId];
79 int nbSize = list_i.size();
80 for (int ne = 0; ne < nbSize; ne++)
81 {
82 int j = list_i[ne];
83 Real r = (pos_i - position[j]).norm();;
84
85 if (r > EPSILON)
86 {
87 Real a_ij = kernWeight(r, smoothingLength);
88 alpha_i += a_ij;
89 }
90 }
91
92 alpha[pId] = alpha_i;
93 }
94
95 template <typename Real>
96 __global__ void SIMPLE_CorrectAlpha
97 (
98 DArray<Real> alpha,
99 Real maxAlpha)
100 {
101 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
102 if (pId >= alpha.size()) return;
103
104 Real alpha_i = alpha[pId];
105 if (alpha_i < maxAlpha)
106 {
107 alpha_i = maxAlpha;
108 }
109 alpha[pId] = alpha_i;
110 }
111
112 template <typename Real, typename Coord>
113 __global__ void SIMPLE_ComputeDiagonalElement
114 (
115 DArray<Real> diaA,
116 DArray<Real> alpha,
117 DArray<Coord> position,
118 DArray<Attribute> attribute,
119 DArrayList<int> neighbors,
120 Real smoothingLength)
121 {
122 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
123 if (pId >= position.size()) return;
124 if (!attribute[pId].isFluid()) return;
125
126 Coord pos_i = position[pId];
127 Real invAlpha_i = 1.0f / alpha[pId];
128 Real A_i = 0.0f;
129
130 List<int>& list_i = neighbors[pId];
131 int nbSize = list_i.size();
132 for (int ne = 0; ne < nbSize; ne++)
133 {
134 int j = list_i[ne];
135 Real r = (pos_i - position[j]).norm();
136
137 if (r > EPSILON && attribute[j].isFluid())
138 {
139 Real wrr_ij = invAlpha_i * kernWRR(r, smoothingLength);
140 A_i += wrr_ij;
141 atomicAdd(&diaA[j], wrr_ij);
142 }
143 }
144
145 atomicAdd(&diaA[pId], A_i);
146 }
147
148 template <typename Real, typename Coord>
149 __global__ void SIMPLE_ComputeDiagonalElement
150 (
151 DArray<Real> AiiFluid,
152 DArray<Real> AiiTotal,
153 DArray<Real> alpha,
154 DArray<Coord> position,
155 DArray<Attribute> attribute,
156 DArrayList<int> neighbors,
157 Real smoothingLength
158 )
159 {
160 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
161 if (pId >= position.size()) return;
162 if (!attribute[pId].isFluid()) return;
163
164 Real invAlpha = 1.0f / alpha[pId];
165
166
167 Real diaA_total = 0.0f;
168 Real diaA_fluid = 0.0f;
169 Coord pos_i = position[pId];
170
171 bool bNearWall = false;
172
173 List<int>& list_i = neighbors[pId];
174 int nbSize = list_i.size();
175 for (int ne = 0; ne < nbSize; ne++)
176 {
177 int j = list_i[ne];
178 //int j = neighbors.getElement(pId, ne);
179 Real r = (pos_i - position[j]).norm();
180
181 Attribute att_j = attribute[j];
182
183
184 if (r > EPSILON)
185 {
186 Real wrr_ij = invAlpha * kernWRR(r, smoothingLength);
187 if (att_j.isFluid())
188 {
189 diaA_total += wrr_ij;
190 diaA_fluid += wrr_ij;
191 atomicAdd(&AiiFluid[j], wrr_ij);
192 atomicAdd(&AiiTotal[j], wrr_ij);
193 }
194 else
195 {
196 diaA_total += 2.0f * wrr_ij;
197 }
198
199 }
200
201 }
202 atomicAdd(&AiiFluid[pId], diaA_fluid);
203 atomicAdd(&AiiTotal[pId], diaA_total);
204
205
206 }
207
208 template <typename Real, typename Coord>
209 __global__ void SIMPLE_DetectSurface
210 (
211 DArray<Real> Aii,
212 DArray<bool> bSurface,
213 DArray<Real> AiiFluid,
214 DArray<Real> AiiTotal,
215 DArray<Coord> position,
216 DArray<Attribute> attribute,
217 DArrayList<int> neighbors,
218 Real smoothingLength,
219 Real maxA
220 )
221 {
222 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
223 if (pId >= position.size()) return;
224 if (!attribute[pId].isFluid()) return;
225
226 Real total_weight = 0.0f;
227 Coord div_i = Coord(0);
228
229 SmoothKernel<Real> kernSmooth;
230
231 Coord pos_i = position[pId];
232
233 List<int>& list_i = neighbors[pId];
234 int nbSize = list_i.size();
235 bool bNearWall = false;
236 for (int ne = 0; ne < nbSize; ne++)
237 {
238
239 int j = list_i[ne];
240 Real r = (pos_i - position[j]).norm();
241
242 if (r > EPSILON && attribute[j].isFluid())
243 {
244 float weight = -kernSmooth.Gradient(r, smoothingLength);
245 total_weight += weight;
246 div_i += (position[j] - pos_i) * (weight / r);
247 }
248
249 if (!attribute[j].isFluid())
250 {
251 bNearWall = true;
252 }
253 }
254
255 total_weight = total_weight < EPSILON ? 1.0f : total_weight;
256 Real absDiv = div_i.norm() / total_weight;
257
258 bool bSurface_i = false;
259 Real diagF_i = AiiFluid[pId];
260 Real diagT_i = AiiTotal[pId];
261 Real aii = diagF_i;
262 Real eps = 0.001f;
263 Real diagS_i = diagT_i - diagF_i;
264 Real threshold = 0.0f;
265 if (bNearWall && diagT_i < maxA * (1.0f - threshold))
266 {
267 bSurface_i = true;
268 aii = maxA - (diagT_i - diagF_i);
269 }
270
271 if (!bNearWall && diagF_i < maxA * (1.0f - threshold))
272 {
273 bSurface_i = true;
274 aii = maxA;
275 }
276 bSurface[pId] = bSurface_i;
277 Aii[pId] = aii;
278 }
279
280
281 template <typename Real, typename Coord>
282 __global__ void SIMPLE_ComputeDivergence
283 (
284 DArray<Real> divergence,
285 DArray<Real> alpha,
286 DArray<Coord> position,
287 DArray<Coord> velocity,
288 DArray<Coord> normals,
289 DArray<Attribute> attribute,
290 DArrayList<int> neighbors,
291 Real separation,
292 Real tangential,
293 Real restDensity,
294 Real smoothingLength,
295 Real dt
296 )
297 {
298 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
299 if (pId >= position.size()) return;
300 if (!attribute[pId].isFluid()) return;
301 Coord pos_i = position[pId];
302 Coord vel_i = velocity[pId];
303
304 Real div_vi = 0.0f;
305
306 Real invAlpha_i = 1.0f / alpha[pId];
307
308 List<int>& list_i = neighbors[pId];
309 int nbSize = list_i.size();
310
311 for (int ne = 0; ne < nbSize; ne++)
312 {
313 int j = list_i[ne];
314
315 Real r = (pos_i - position[j]).norm();
316
317 if (r > EPSILON && attribute[j].isFluid())
318 {
319 Real wr_ij = kernWR(r, smoothingLength);
320 Coord g = -invAlpha_i * (pos_i - position[j]) * wr_ij * (1.0f / r);
321
322 if (attribute[j].isFluid())
323 {
324 Real div_ij = 0.5f * (vel_i - velocity[j]).dot(g) * restDensity / dt; //dv_ij = 1 / alpha_i * (v_i-v_j).*(x_i-x_j) / r * (w / r);
325 atomicAdd(&divergence[pId], div_ij);
326 atomicAdd(&divergence[j], div_ij);
327 }
328 else
329 {
330
331 Coord normal_j = normals[j];
332
333 Coord dVel = vel_i - velocity[j];
334 Real magNVel = dVel.dot(normal_j);
335 Coord nVel = magNVel * normal_j;
336 Coord tVel = dVel - nVel;
337
338
339 if (magNVel < -EPSILON)
340 {
341 Real div_ij = g.dot(2.0f * (nVel + tangential * tVel)) * restDensity / dt;
342 // printf("Boundary div: %f \n", div_ij);
343 atomicAdd(&divergence[pId], div_ij);
344 }
345 else
346 {
347 Real div_ij = g.dot(2.0f * (separation * nVel + tangential * tVel)) * restDensity / dt;
348 atomicAdd(&divergence[pId], div_ij);
349 }
350
351 }
352 }
353 }
354 }
355
356 template <typename Real, typename Coord>
357 __global__ void SIMPLE_CompensateSource
358 (
359 DArray<Real> divergence,
360 DArray<Real> density,
361 DArray<Attribute> attribute,
362 DArray<Coord> position,
363 Real restDensity,
364 Real dt
365 )
366 {
367 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
368 if (pId >= density.size()) return;
369 if (!attribute[pId].isFluid()) return;
370
371 Coord pos_i = position[pId];
372 if (density[pId] > restDensity)
373 {
374 Real ratio = (density[pId] - restDensity) / restDensity;
375 atomicAdd(&divergence[pId], 100000.0f * ratio / dt);
376 }
377 }
378
379
380 // compute Ax;
381 template <typename Real, typename Coord>
382 __global__ void SIMPLE_ComputeAx
383 (
384 DArray<Real> residual,
385 DArray<Real> pressure,
386 DArray<Real> aiiSymArr,
387 DArray<Real> alpha,
388 DArray<Coord> position,
389 DArray<Attribute> attribute,
390 DArrayList<int> neighbor,
391 Real smoothingLength
392 )
393 {
394 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
395 if (pId >= position.size()) return;
396 if (!attribute[pId].isFluid()) return;
397
398 Coord pos_i = position[pId];
399 Real invAlpha_i = 1.0f / alpha[pId];
400
401 atomicAdd(&residual[pId], aiiSymArr[pId] * pressure[pId]);
402 Real con1 = 1.0f;// PARAMS.mass / PARAMS.restDensity / PARAMS.restDensity;
403
404 //int nbSize = neighbor.getNeighborSize(pId);
405 List<int>& list_i = neighbor[pId];
406 int nbSize = list_i.size();
407 for (int ne = 0; ne < nbSize; ne++)
408 {
409 int j = list_i[ne];
410 //int j = neighbor.getElement(pId, ne);
411
412 Real r = (pos_i - position[j]).norm();
413
414 if (r > EPSILON && attribute[j].isFluid())
415 {
416 Real wrr_ij = kernWRR(r, smoothingLength);
417 Real a_ij = -invAlpha_i * wrr_ij;
418 // residual += con1*a_ij*preArr[j];
419 atomicAdd(&residual[pId], con1 * a_ij * pressure[j]);
420 atomicAdd(&residual[j], con1 * a_ij * pressure[pId]);
421 }
422 }
423 }
424
425
426 template <typename Real, typename Coord>
427 __global__ void SIMPLE_P_dv(
428 DArray<Coord> P_dv,
429 DArray<Real> pressure,
430 DArray<Real> alpha,
431 DArray<bool> bSurface,
432 DArray<Coord> position,
433 DArray<Coord> velocity,
434 DArray<Coord> normal,
435 DArray<Attribute> attribute,
436 DArrayList<int> neighbor,
437 Real restDensity,
438 Real airPressure,
439 Real sliding,
440 Real separation,
441 Real smoothingLength,
442 Real dt)
443 {
444 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
445 if (pId >= position.size()) return;
446
447 if (attribute[pId].isFluid())
448 {
449 Coord pos_i = position[pId];
450 Real p_i = pressure[pId];
451
452 //int nbSize = neighbor.getNeighborSize(pId);
453
454 Real total_weight = 0.0f;
455
456 Real ceo = 1.6f;
457
458 Real invAlpha = 1.0f / alpha[pId];
459 Coord vel_i = velocity[pId];
460 Coord dv_i(0.0f);
461 Real scale = 1.0f * dt / restDensity;
462 Real acuP = 0.0f;
463 total_weight = 0.0f;
464
465
466 List<int>& list_i = neighbor[pId];
467 int nbSize = list_i.size();
468 for (int ne = 0; ne < nbSize; ne++)
469 {
470 //int j = neighbor.getElement(pId, ne);
471 int j = list_i[ne];
472 Real r = (pos_i - position[j]).norm();
473
474 Attribute att_j = attribute[j];
475 if (r > EPSILON)
476 {
477 Real weight = -invAlpha * kernWR(r, smoothingLength);
478 Coord dnij = (pos_i - position[j]) * (1.0f / r);
479 Coord corrected = dnij;
480 if (corrected.norm() > EPSILON)
481 {
482 corrected = corrected.normalize();
483 }
484 corrected = -scale * weight * corrected;
485
486 Coord dvij = (pressure[j] - pressure[pId]) * corrected;
487 Coord dvjj = (pressure[j] + airPressure) * corrected;
488 Coord dvij_sym = 0.5f * (pressure[pId] + pressure[j]) * corrected;
489
490
491 if (att_j.isFluid())
492 {
493 if (bSurface[pId])
494 {
495 dv_i += dvjj;
496 }
497 else
498 {
499 dv_i += dvij;
500 }
501
502 if (bSurface[j])
503 {
504 Coord dvii = -(pressure[pId] + airPressure) * corrected;
505 atomicAdd(&P_dv[j][0], ceo * dvii[0]);
506 atomicAdd(&P_dv[j][1], ceo * dvii[1]);
507 atomicAdd(&P_dv[j][2], ceo * dvii[2]);
508 }
509 else
510 {
511
512 atomicAdd(&P_dv[j][0], ceo * dvij[0]);
513 atomicAdd(&P_dv[j][1], ceo * dvij[1]);
514 atomicAdd(&P_dv[j][2], ceo * dvij[2]);
515 }
516 }
517 else
518 {
519 Coord dvii = 2.0f * (pressure[pId]) * corrected;
520 if (bSurface[pId])
521 {
522 dv_i += dvii;
523 }
524
525 float weight = 2.0f * invAlpha * kernWeight(r, smoothingLength);
526 Coord nij = (pos_i - position[j]);
527 if (nij.norm() > EPSILON)
528 {
529 nij = nij.normalize();
530 }
531 else
532 nij = Coord(1.0f, 0.0f, 0.0f);
533
534 Coord normal_j = normal[j];
535 Coord dVel = velocity[j] - vel_i;
536 Real magNVel = dVel.dot(normal_j);
537 Coord nVel = magNVel * normal_j;
538 Coord tVel = dVel - nVel;
539 if (magNVel > EPSILON)
540 {
541 dv_i += weight * nij.dot(nVel + sliding * tVel) * nij;
542 }
543 else
544 {
545 dv_i += weight * nij.dot(separation * nVel + sliding * tVel) * nij;
546 }
547 }
548 }
549 }
550 dv_i *= ceo;
551 atomicAdd(&P_dv[pId][0], dv_i[0]);
552 atomicAdd(&P_dv[pId][1], dv_i[1]);
553 atomicAdd(&P_dv[pId][2], dv_i[2]);
554 }
555 }
556
557
558
559 template <typename Real>
560 __global__ void UpdatePressure(
561 DArray<Real> totalPressure,
562 DArray<Real> deltaPressure)
563 {
564 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
565 if (pId >= totalPressure.size()) return;
566
567 totalPressure[pId] += deltaPressure[pId];
568 }
569
570
571 template <typename Coord>
572 __global__ void SIMPLE_VelUpdate(
573 DArray<Coord> velocity,
574 DArray<Coord> vel_old,
575 DArray<Coord> P_dv
576 )
577 {
578 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
579 if (pId >= velocity.size()) return;
580 velocity[pId] = vel_old[pId] + P_dv[pId];
581 }
582
583
584 template <typename Real, typename Coord, typename Matrix>
585 __global__ void SIMPLE_VisComput
586 (
587 DArray<Coord> velNew,
588 DArray<Coord> velBuf,
589 DArray<Coord> velOld,
590 DArray<Coord> velDp,
591 DArray<Coord> position,
592 DArray<Real> alpha,
593 DArray<Attribute> attribute,
594 DArrayList<int> neighbor,
595 Real rest_density,
596 Real h,
597 Real dt,
598 DArray<Real> vis
599 )
600 {
601 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
602 if (pId >= position.size()) return;
603 if (!attribute[pId].isFluid()) return;
604
605 Real tempValue = 0;
606 Coord Avj(0);
607 Matrix Mii(0);
608 Real invAlpha_i = 1.0f / alpha[pId];
609 //int nbSize = neighbor.getNeighborSize(pId);
610 List<int>& list_i = neighbor[pId];
611 int nbSize = list_i.size();
612
613 for (int ne = 0; ne < nbSize; ne++)
614 {
615 int j = list_i[ne];
616 //int j = neighbor.getElement(pId, ne);
617 Real r = (position[pId] - position[j]).norm();
618 if (r > EPSILON)
619 {
620 Real wrr_ij = kernWRR(r, h);
621 Real ai = 1.0f / alpha[pId];
622 Real aj = 1.0f / alpha[j];
623 Coord nij = (position[j] - position[pId]) / r;
624 tempValue = 0.25 * dt * (vis[pId] + vis[j]) * (ai + aj) * wrr_ij / rest_density;
625 Avj += tempValue * velBuf[j].dot(nij) * nij;
626 Matrix Mij(0);
627 Mij(0, 0) += nij[0] * nij[0]; Mij(0, 1) += nij[0] * nij[1]; Mij(0, 2) += nij[0] * nij[2];
628 Mij(1, 0) += nij[1] * nij[0]; Mij(1, 1) += nij[1] * nij[1]; Mij(1, 2) += nij[1] * nij[2];
629 Mij(2, 0) += nij[2] * nij[0]; Mij(2, 1) += nij[2] * nij[1]; Mij(2, 2) += nij[2] * nij[2];
630 Mii += Mij * tempValue;
631 }
632 }
633 Mii += Matrix::identityMatrix();
634 velNew[pId] = Mii.inverse() * (velOld[pId] + velDp[pId] + Avj);
635 }
636
637 template <typename Real, typename Coord>
638 __global__ void SIMPLE_Vis_AxComput
639 (
640 DArray<Real> v_y,
641 DArray<Coord> velBuf,
642 DArray<Coord> position,
643 DArray<Real> alpha,
644 DArray<Attribute> attribute,
645 DArrayList<int> neighbor,
646 Real rest_density,
647 Real h,
648 Real dt,
649 DArray<Real> vis
650 )
651 {
652 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
653 if (pId >= position.size()) return;
654 if (!attribute[pId].isFluid()) return;
655
656 Real tempValue = 0;
657 Coord Avi(0);
658 Real invAlpha_i = 1.0f / alpha[pId];
659 //int nbSize = neighbor.getNeighborSize(pId);
660 List<int>& list_i = neighbor[pId];
661 int nbSize = list_i.size();
662 for (int ne = 0; ne < nbSize; ne++)
663 {
664 int j = list_i[ne];
665 //int j = neighbor.getElement(pId, ne);
666 Real r = (position[pId] - position[j]).norm();
667 if (r > EPSILON)
668 {
669 Real wrr_ij = kernWRR(r, h);
670 Real ai = 1.0f / alpha[pId];
671 Real aj = 1.0f / alpha[j];
672 Coord nij = (position[j] - position[pId]) / r;
673 Coord vij = velBuf[pId] - velBuf[j];
674 tempValue = 2.5 * dt * (vis[pId] + vis[j]) * (ai + aj) * wrr_ij / rest_density;
675 Avi += tempValue * vij.dot(nij) * nij;
676 }
677 }
678 Avi += velBuf[pId];
679 v_y[3 * pId] = Avi[0];
680 v_y[3 * pId + 1] = Avi[1];
681 v_y[3 * pId + 2] = Avi[2];
682
683
684 }
685
686
687 template <typename Real, typename Coord>
688 __global__ void SIMPLE_Vis_r_Comput
689 (
690 DArray<Real> v_r,
691 DArray<Real> v_y,
692 DArray<Coord> vel_old,
693 DArray<Coord> Dvel,
694 DArray<Coord> position,
695 DArray<Attribute> attribute
696 )
697 {
698 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
699 if (pId >= position.size()) return;
700 if (!attribute[pId].isFluid()) return;
701 Coord temp_vel = Dvel[pId] + vel_old[pId];
702 v_r[3 * pId] = temp_vel[0] - v_y[3 * pId];
703 v_r[3 * pId + 1] = temp_vel[1] - v_y[3 * pId + 1];
704 v_r[3 * pId + 2] = temp_vel[2] - v_y[3 * pId + 2];
705 }
706
707 template <typename Real, typename Coord>
708 __global__ void SIMPLE_Vis_pToVector
709 (
710 DArray<Real> v_p,
711 DArray<Coord> pv,
712 DArray<Attribute> attribute
713 )
714 {
715 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
716 if (pId >= attribute.size()) return;
717 if (!attribute[pId].isFluid()) return;
718 pv[pId][0] = v_p[3 * pId];
719 pv[pId][1] = v_p[3 * pId + 1];
720 pv[pId][2] = v_p[3 * pId + 2];
721 }
722
723
724 template <typename Real, typename Coord>
725 __global__ void SIMPLE_Vis_CoordToReal
726 (
727 DArray<Real> veloReal,
728 DArray<Coord> vel,
729 DArray<Attribute> attribute
730 )
731 {
732 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
733 if (pId >= attribute.size()) return;
734 if (!attribute[pId].isFluid()) return;
735 veloReal[3 * pId] = vel[pId][0];
736 veloReal[3 * pId + 1] = vel[pId][1];
737 veloReal[3 * pId + 2] = vel[pId][2];
738 }
739
740 template <typename Real, typename Coord>
741 __global__ void SIMPLE_Vis_RealToVeloctiy
742 (
743 DArray<Coord> vel,
744 DArray<Real> veloReal,
745 DArray<Attribute> attribute
746 )
747 {
748 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
749 if (pId >= attribute.size()) return;
750 if (!attribute[pId].isFluid()) return;
751 vel[pId][0] = veloReal[3 * pId];
752 vel[pId][1] = veloReal[3 * pId + 1];
753 vel[pId][2] = veloReal[3 * pId + 2];
754 }
755
756 template <typename Real, typename Coord>
757 __global__ void SIMPLE_CrossVis
758 (
759 DArray<Real> crossVis,
760 DArray<Coord> velBuf,
761 DArray<Coord> position,
762 DArray<Real> alpha,
763 DArray<Attribute> attribute,
764 DArrayList<int> neighbor,
765 Real smoothingLength,
766 Real visTop,
767 Real visFloor,
768 Real K,
769 Real N
770 )
771 {
772 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
773 if (pId >= position.size()) return;
774 if (!attribute[pId].isFluid()) return;
775
776 Coord pos_i = position[pId];
777 Coord vel_i = velBuf[pId];
778 Real invAlpha_i = 1.0f / alpha[pId];
779
780 Coord dv(0);
781 Coord vij(0);
782
783 List<int>& list_i = neighbor[pId];
784 int nbSize = list_i.size();
785 //int nbSize = neighbor.getNeighborSize(pId);
786 for (int ne = 0; ne < nbSize; ne++)
787 {
788 //int j = neighbor.getElement(pId, ne);
789 int j = list_i[ne];
790 Real r = (position[pId] - position[j]).norm();
791 if (r > EPSILON && attribute[j].isFluid())
792 {
793 Real wr_ij = kernWR(r, smoothingLength);
794 Coord g = -invAlpha_i * (pos_i - position[j]) * wr_ij * (1.0f / r);
795 vij = 0.5 * (velBuf[pId] - velBuf[j]);
796 dv[0] += vij[0] * g[0];
797 dv[1] += vij[1] * g[1];
798 dv[2] += vij[2] * g[2];
799 }
800
801 }
802 Real Norm = VB_FNorm(dv[0], dv[1], dv[2]);
803 crossVis[pId] = VB_Viscosity(visTop, visFloor, Norm, K, N);
804 if (pId == 100) printf("viscosity : %f. \r\n", crossVis[pId]);
805 }
806
807
808 __global__ void SIMPLE_AttributeInit
809 (
810 DArray<Attribute> atts
811 )
812 {
813 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
814 if (pId >= atts.size()) return;
815
816 atts[pId].setFluid();
817 atts[pId].setDynamic();
818 }
819
820
821 template<typename TDataType>
822 SimpleVelocityConstraint<TDataType>::SimpleVelocityConstraint()
823 : ConstraintModule()
824 , m_airPressure(Real(0))
825 , m_reduce(NULL)
826 , m_arithmetic(NULL)
827 {
828 this->inSmoothingLength()->setValue(Real(0.0125));
829 this->varRestDensity()->setValue(Real(1000));
830
831 m_densitySum = std::make_shared<SummationDensity<TDataType>>();
832 this->varRestDensity()->connect(m_densitySum->varRestDensity());
833 this->inSmoothingLength()->connect(m_densitySum->inSmoothingLength());
834 this->inSamplingDistance()->connect(m_densitySum->inSamplingDistance());
835
836 this->inPosition()->connect(m_densitySum->inPosition());
837 this->inNeighborIds()->connect(m_densitySum->inNeighborIds());
838 SIMPLE_IterNum = 5;
839
840 this->inAttribute()->tagOptional(true);
841 this->inNormal()->tagOptional(true);
842 };
843
844 template<typename TDataType>
845 SimpleVelocityConstraint<TDataType>::~SimpleVelocityConstraint()
846 {
847 m_alpha.clear();
848 m_Aii.clear();
849 m_AiiFluid.clear();
850 m_AiiTotal.clear();
851 m_pressure.clear();
852 m_divergence.clear();
853 m_bSurface.clear();
854
855 m_y.clear();
856 m_r.clear();
857 m_p.clear();
858
859 m_pressure.clear();
860
861 if (m_reduce)
862 {
863 delete m_reduce;
864 }
865 if (m_arithmetic)
866 {
867 delete m_arithmetic;
868 }
869 };
870
871
872 template<typename TDataType>
873 void SimpleVelocityConstraint<TDataType>::constrain()
874 {
875
876 if ((init_flag == false) || (this->inTimeStep()->getValue() == 0))
877 {
878 initialize();
879 }
880
881 if (this->inPosition()->size() != m_viscosity.size())
882 {
883 resizeVector();
884 }
885
886 cuSynchronize();
887
888 auto& m_position = this->inPosition()->getData();
889 auto& m_velocity = this->inVelocity()->getData();
890 auto& m_neighborhood = this->inNeighborIds()->getData();
891 auto& m_attribute = this->inAttribute()->getData();
892 auto& m_normal = this->inNormal()->getData();
893 auto m_smoothingLength = this->inSmoothingLength()->getValue();
894 auto m_restDensity = this->varRestDensity()->getValue();
895
896 //Real dt = getParent()->getDt();
897 Real dt = this->inTimeStep()->getData();
898
899 int num = this->inPosition()->size();
900 uint pDims = cudaGridSize(this->inPosition()->size(), BLOCK_SIZE);
901
902
903 m_alpha.reset();
904 //compute alpha_i = sigma w_j and A_i = sigma w_ij / r_ij / r_ij
905 SIMPLE_ComputeAlpha << <pDims, BLOCK_SIZE >> > (
906 m_alpha,
907 m_position,
908 m_attribute,
909 m_neighborhood,
910 m_smoothingLength);
911
912 SIMPLE_CorrectAlpha << <pDims, BLOCK_SIZE >> > (
913 m_alpha,
914 m_maxAlpha);
915
916
917 //compute the diagonal elements of the coefficient matrix
918 m_AiiFluid.reset();
919 m_AiiTotal.reset();
920
921 SIMPLE_ComputeDiagonalElement << <pDims, BLOCK_SIZE >> > (
922 m_AiiFluid,
923 m_AiiTotal,
924 m_alpha,
925 m_position,
926 m_attribute,
927 m_neighborhood,
928 m_smoothingLength);
929
930 m_bSurface.reset();
931 m_Aii.reset();
932 SIMPLE_DetectSurface << <pDims, BLOCK_SIZE >> > (
933 m_Aii,
934 m_bSurface,
935 m_AiiFluid,
936 m_AiiTotal,
937 m_position,
938 m_attribute,
939 m_neighborhood,
940 m_smoothingLength,
941 m_maxA);
942
943
944 m_densitySum->compute();
945 m_densitySum->outDensity()->getData();
946
947 //IF Cross Model is active, viscous coefficients should be computed.
948 if (IsCrossReady)
949 {
950
951 SIMPLE_CrossVis << <pDims, BLOCK_SIZE >> > (
952 m_viscosity,
953 m_velocity,
954 m_position,
955 m_alpha,
956 m_attribute,
957 m_neighborhood,
958 m_smoothingLength,
959 CrossVisCeil,
960 CrossVisFloor,
961 Cross_K,
962 Cross_N
963 );
964 }
965
966
967 int totalIter = 0;
968
969 m_pressure.reset();
970
971 velOld.assign(m_velocity);
972
973 Real Old_temp = 0;
974
975 if (this->varSimpleIterationEnable()->getValue() == false)
976 {
977 SIMPLE_IterNum = 1;
978 }
979
980 //SIMPLE Algorithm / Outer Iterations
981 while (totalIter < SIMPLE_IterNum)
982 {
983 printf("Iteration : %d *****", totalIter);
984 totalIter++;
985
986 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
987 //Incopressibility Solver -- Begin ¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý
988 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
989 m_divergence.reset();
990 SIMPLE_ComputeDivergence << <pDims, BLOCK_SIZE >> > (
991 m_divergence,
992 m_alpha,
993 m_position,
994 m_velocity,
995 m_normal,
996 m_attribute,
997 m_neighborhood,
998 m_separation,
999 m_tangential,
1000 m_restDensity,
1001 m_smoothingLength,
1002 dt);
1003
1004 SIMPLE_CompensateSource << <pDims, BLOCK_SIZE >> > (
1005 m_divergence,
1006 m_densitySum->outDensity()->getData(),
1007 m_attribute,
1008 m_position,
1009 m_restDensity,
1010 dt);
1011
1012 m_deltaPressure.reset();
1013 m_y.reset();
1014 SIMPLE_ComputeAx << <pDims, BLOCK_SIZE >> > (
1015 m_y,
1016 m_deltaPressure,
1017 m_Aii,
1018 m_alpha,
1019 m_position,
1020 m_attribute,
1021 m_neighborhood,
1022 m_smoothingLength);
1023
1024 m_r.reset();
1025 Function2Pt::subtract(m_r, m_divergence, m_y);
1026
1027 m_p.assign(m_r);
1028 Real rr = m_arithmetic->Dot(m_r, m_r);
1029 Real err = sqrt(rr / m_r.size());
1030
1031 Real initErr = err;
1032
1033 int itor = 0;
1034 while (itor < 1000 && err / initErr > 0.00001f)
1035 {
1036 m_y.reset();
1037 SIMPLE_ComputeAx << <pDims, BLOCK_SIZE >> > (
1038 m_y,
1039 m_p,
1040 m_Aii,
1041 m_alpha,
1042 m_position,
1043 m_attribute,
1044 m_neighborhood,
1045 m_smoothingLength);
1046
1047 float alpha = rr / m_arithmetic->Dot(m_p, m_y);
1048 Function2Pt::saxpy(m_deltaPressure, m_p, m_deltaPressure, alpha);
1049 Function2Pt::saxpy(m_r, m_y, m_r, -alpha);
1050
1051 Real rr_old = rr;
1052
1053 rr = m_arithmetic->Dot(m_r, m_r);
1054
1055 Real beta = rr / rr_old;
1056 Function2Pt::saxpy(m_p, m_p, m_r, beta);
1057
1058 err = sqrt(rr / m_r.size());
1059 itor++;
1060 }
1061
1062 UpdatePressure << <pDims, BLOCK_SIZE >> > (
1063 m_pressure,
1064 m_deltaPressure);
1065
1066 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1067 //Incopressibility Solver -- End ¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü
1068 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1069
1070
1071
1072 P_dv.reset();
1073 SIMPLE_P_dv << <pDims, BLOCK_SIZE >> > (
1074 P_dv,
1075 m_pressure,
1076 m_alpha,
1077 m_bSurface,
1078 m_position,
1079 m_velocity,
1080 m_normal,
1081 m_attribute,
1082 m_neighborhood,
1083 m_restDensity,
1084 m_airPressure,
1085 m_tangential,
1086 m_separation,
1087 m_smoothingLength,
1088 dt);
1089
1090 SIMPLE_VelUpdate << <pDims, BLOCK_SIZE >> > (
1091 m_velocity,
1092 velOld,
1093 P_dv
1094 );
1095
1096
1097 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1098 //Viscosity Solver -- Begin ¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý¡ý
1099 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1100 v_y.reset();
1101 SIMPLE_Vis_AxComput<Real, Coord> << <pDims, BLOCK_SIZE >> > (
1102 v_y,
1103 m_velocity,
1104 m_position,
1105 m_alpha,
1106 m_attribute,
1107 m_neighborhood,
1108 m_restDensity,
1109 m_smoothingLength,
1110 dt,
1111 m_viscosity
1112 );
1113
1114 v_r.reset();
1115 SIMPLE_Vis_r_Comput << <pDims, BLOCK_SIZE >> > (
1116 v_r,
1117 v_y,
1118 velOld,
1119 P_dv,
1120 m_position,
1121 m_attribute
1122 );
1123 v_p.assign(v_r);
1124
1125
1126 Real Vrr = m_arithmetic_v->Dot(v_r, v_r);
1127 Real Verr = sqrt(Vrr / v_r.size());
1128 int VisItor = 0;
1129 initErr = Verr;
1130
1131 while (VisItor < 1000 && Verr / initErr > 0.00001f)
1132 {
1133 VisItor++;
1134 //The type of "v_p" should convert to DArray<Coord>
1135 SIMPLE_Vis_pToVector<Real, Coord> << <pDims, BLOCK_SIZE >> > (
1136 v_p,
1137 v_pv,
1138 m_attribute
1139 );
1140
1141 v_y.reset();
1142 SIMPLE_Vis_AxComput<Real, Coord> << <pDims, BLOCK_SIZE >> > (
1143 v_y,
1144 v_pv,
1145 //m_velocity.getValue(),
1146 m_position,
1147 m_alpha,
1148 m_attribute,
1149 m_neighborhood,
1150 m_restDensity,
1151 m_smoothingLength,
1152 dt,
1153 m_viscosity
1154 );
1155
1156 float alpha = Vrr / m_arithmetic_v->Dot(v_p, v_y);
1157
1158 //The type of "velocity" should convert to DArray<Real>
1159 SIMPLE_Vis_CoordToReal<Real, Coord> << <pDims, BLOCK_SIZE >> > (
1160 m_VelocityReal,
1161 m_velocity,
1162 m_attribute
1163 );
1164
1165 Function2Pt::saxpy(m_VelocityReal, v_p, m_VelocityReal, alpha);
1166
1167 SIMPLE_Vis_RealToVeloctiy<Real, Coord> << <pDims, BLOCK_SIZE >> > (
1168 m_velocity,
1169 m_VelocityReal,
1170 m_attribute
1171 );
1172
1173 Function2Pt::saxpy(v_r, v_y, v_r, -alpha);
1174 Real Vrr_old = Vrr;
1175
1176 Vrr = m_arithmetic_v->Dot(v_r, v_r);
1177 Real beta = Vrr / Vrr_old;
1178 Function2Pt::saxpy(v_p, v_p, v_r, beta);
1179
1180 Verr = sqrt(Vrr / v_r.size());
1181
1182 }
1183 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1184 //Viscosity Solver -- End ¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü¡ü
1185 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1186
1187 Real ppp = m_arithmetic->Dot(m_pressure, m_pressure);
1188 ppp = sqrt(ppp / m_pressure.size());
1189 Real testP = (ppp - Old_temp) ;
1190 printf("%f \r\n", testP);
1191 Old_temp = ppp;
1192 }
1193 };
1194
1195 template<typename TDataType>
1196 void SimpleVelocityConstraint<TDataType>::initialAttributes() {
1197
1198 if (this->inAttribute()->isEmpty() || this->inAttribute()->size() != this->inPosition()->size())
1199 {
1200 this->inAttribute()->allocate();
1201 this->inAttribute()->resize(this->inPosition()->size());
1202 this->inNormal()->allocate();
1203 this->inNormal()->resize(this->inPosition()->size());
1204 this->inNormal()->getData().reset();
1205 cuExecute(this->inPosition()->size(), SIMPLE_AttributeInit, this->inAttribute()->getData());
1206 }
1207 }
1208
1209
1210 template<typename TDataType>
1211 bool SimpleVelocityConstraint<TDataType>::resizeVector()
1212 {
1213
1214 int num = this->inPosition()->size();
1215 m_alpha.resize(num);
1216 m_Aii.resize(num);
1217 m_AiiFluid.resize(num);
1218 m_AiiTotal.resize(num);
1219 m_pressure.resize(num);
1220 m_divergence.resize(num);
1221 m_bSurface.resize(num);
1222
1223
1224 m_y.resize(num);
1225 m_r.resize(num);
1226 m_p.resize(num);
1227
1228 v_y.resize(3 * num);
1229 v_r.resize(3 * num);
1230 v_p.resize(3 * num);
1231 v_pv.resize(num);
1232 m_VelocityReal.resize(3 * num);
1233
1234 m_pressure.resize(num);
1235
1236 P_dv.resize(num);
1237 velOld.resize(num);
1238 velBuf.resize(num);
1239 m_viscosity.resize(num);
1240 m_deltaPressure.resize(num);
1241 m_pressBuf.resize(num);
1242 m_crossViscosity.resize(num);
1243
1244
1245 m_reduce = Reduction<float>::Create(num);
1246 m_arithmetic = Arithmetic<float>::Create(num);
1247 m_arithmetic_v = Arithmetic<float>::Create(3 * num);
1248
1249 visValueSet();
1250
1251 initialAttributes();
1252
1253 return true;
1254 }
1255
1256
1257
1258 template<typename TDataType>
1259 bool SimpleVelocityConstraint<TDataType>::initialize()
1260 {
1261 cuSynchronize();
1262 init_flag = true;
1263
1264 initialAttributes();
1265
1266 auto& m_position = this->inPosition()->getData();
1267 auto& m_velocity = this->inVelocity()->getData();
1268 auto& m_neighborhood = this->inNeighborIds()->getData();
1269 auto& m_attribute = this->inAttribute()->getData();
1270 auto m_smoothingLength = this->inSmoothingLength()->getValue();
1271 auto m_restDensity = this->varRestDensity()->getValue();
1272 auto& m_normal = this->inNormal()->getData();
1273
1274 resizeVector();
1275
1276
1277 uint pDims = cudaGridSize(this->inPosition()->size(), BLOCK_SIZE);
1278
1279 m_alpha.reset();
1280
1281 SIMPLE_ComputeAlpha << <pDims, BLOCK_SIZE >> > (
1282 m_alpha,
1283 m_position,
1284 m_attribute,
1285 m_neighborhood,
1286 m_smoothingLength);
1287
1288 m_maxAlpha = m_reduce->maximum(m_alpha.begin(), m_alpha.size());
1289
1290 SIMPLE_CorrectAlpha << <pDims, BLOCK_SIZE >> > (
1291 m_alpha,
1292 m_maxAlpha);
1293
1294 m_AiiFluid.reset();
1295 SIMPLE_ComputeDiagonalElement << <pDims, BLOCK_SIZE >> > (
1296 m_AiiFluid,
1297 m_alpha,
1298 m_position,
1299 m_attribute,
1300 m_neighborhood,
1301 m_smoothingLength);
1302
1303 m_maxA = m_reduce->maximum(m_AiiFluid.begin(), m_AiiFluid.size());
1304
1305
1306
1307 std::cout << "Max alpha: " << m_maxAlpha << std::endl;
1308 std::cout << "Max A: " << m_maxA << std::endl;
1309
1310 return true;
1311 };
1312
1313
1314 DEFINE_CLASS(SimpleVelocityConstraint);
1315
1316}
Generated by doxygen 1.13.2