PeriDyno 1.0.0
Loading...
Searching...
No Matches
ApproximateImplicitViscosity.cu
Go to the documentation of this file.
1#include "ApproximateImplicitViscosity.h"
2//#include <string>
3#include "Algorithm/Function2Pt.h"
4
5
6namespace dyno
7{
8
9 template<typename TDataType>
10 ApproximateImplicitViscosity<TDataType>::ApproximateImplicitViscosity()
11 :ConstraintModule()
12 {
13 this->inAttribute()->tagOptional(true);
14 }
15
16
17 //F-norms
18 template<typename Real>
19 __device__ Real VB_FNorm(const Real a, const Real b, const Real c)
20 {
21 Real p = a + b + c;
22 return sqrt(p * p / 2);
23 }
24
25 //CrossModel Viscosity Coefficient
26 template<typename Real>
27 __device__ Real VB_Viscosity(const Real Viscosity_h, const Real Viscosity_l, const Real StrainRate, const Real CrossModel_K, const Real CrossModel_n)
28 {
29 Real p = CrossModel_K * StrainRate;
30 p = pow(p, CrossModel_n);
31 return Viscosity_h + (Viscosity_l - Viscosity_h) / (1 + p);
32 }
33
34 __device__ inline float kernWeight(const float r, const float h)
35 {
36 const float q = r / h;
37 if (q > 1.0f) return 0.0f;
38 else {
39 return (1.0 - pow(q, 4.0f));
40 }
41 }
42
43 __device__ inline float kernWR(const float r, const float h)
44 {
45 float w = kernWeight(r, h);
46 const float q = r / h;
47 if (q < 0.4f)
48 {
49 return w / (0.4f*h);
50 }
51 return w / r;
52 }
53
54 __device__ inline float kernWRR(const float r, const float h)
55 {
56 float w = kernWeight(r, h);
57 const float q = r / h;
58 if (q < 0.4f)
59 {
60 return w / (0.16f*h*h);
61 }
62 return w / r / r;
63 }
64
65 template <typename Real, typename Coord>
66 __global__ void VC_ComputeAlpha
67 (
68 DArray<Real> alpha,
69 DArray<Coord> position,
70 DArray<Attribute> attribute,
71 DArrayList<int> neighbors,
72 Real smoothingLength
73 )
74 {
75 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
76 if (pId >= position.size()) return;
77 if (!attribute[pId].isDynamic()) return;
78
79 Coord pos_i = position[pId];
80 Real alpha_i = 0.0f;
81
82 List<int>& list_i = neighbors[pId];
83 int nbSize = list_i.size();
84
85 for (int ne = 0; ne < nbSize; ne++)
86 {
87 int j = list_i[ne];
88 Real r = (pos_i - position[j]).norm();;
89
90 if (r > EPSILON)
91 {
92 Real a_ij = kernWeight(r, smoothingLength);
93 alpha_i += a_ij;
94 }
95 }
96 alpha[pId] = alpha_i;
97 }
98
99 template <typename Real>
100 __global__ void VC_CorrectAlpha
101 (
102 DArray<Real> alpha,
103 Real maxAlpha)
104 {
105 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
106 if (pId >= alpha.size()) return;
107
108 Real alpha_i = alpha[pId];
109 if (alpha_i < maxAlpha)
110 {
111 alpha_i = maxAlpha;
112 }
113 alpha[pId] = alpha_i;
114 }
115
116 //Jacobi Mehthod
117 template <typename Real, typename Coord, typename Matrix>
118 __global__ void VC_VisComput
119 (
120 DArray<Coord> velNew,
121 DArray<Coord> velBuf,
122 DArray<Coord> velOld,
123 DArray<Coord> velDp,
124 DArray<Coord> position,
125 DArray<Real> alpha,
126 DArray<Attribute> attribute,
127 DArrayList<int> neighbor,
128 Real rest_density,
129 Real h,
130 Real dt,
131 DArray<Real> vis
132 )
133 {
134 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
135 if (pId >= position.size()) return;
136 if (!attribute[pId].isDynamic()) return;
137
138 Real tempValue = 0;
139 Coord Avj(0);
140 Matrix Mii(0);
141 Real invAlpha_i = 1.0f / alpha[pId];
142
143 List<int>& list_i = neighbor[pId];
144 int nbSize = list_i.size();
145 for (int ne = 0; ne < nbSize; ne++)
146 {
147 int j = list_i[ne];
148 Real r = (position[pId] - position[j]).norm();
149 if (r > EPSILON)
150 {
151 Real wrr_ij = kernWRR(r, h);
152 Real ai = 1.0f / alpha[pId];
153 Real aj = 1.0f / alpha[j];
154 Coord nij = (position[j] - position[pId]) / r;
155 tempValue = 0.25 * dt * (vis[pId] + vis[j]) * (ai + aj) * wrr_ij / rest_density;
156 Avj += tempValue * velBuf[j].dot(nij) * nij;
157 Matrix Mij(0);
158 Mij(0, 0) += nij[0] * nij[0]; Mij(0, 1) += nij[0] * nij[1]; Mij(0, 2) += nij[0] * nij[2];
159 Mij(1, 0) += nij[1] * nij[0]; Mij(1, 1) += nij[1] * nij[1]; Mij(1, 2) += nij[1] * nij[2];
160 Mij(2, 0) += nij[2] * nij[0]; Mij(2, 1) += nij[2] * nij[1]; Mij(2, 2) += nij[2] * nij[2];
161 Mii += Mij*tempValue;
162 }
163 }
164 Mii += Matrix::identityMatrix();
165 velNew[pId] = Mii.inverse()*(velOld[pId] + velDp[pId] + Avj);
166 }
167
168
169 //Conjugate Gradient Method.
170 template <typename Real, typename Coord>
171 __global__ void VC_Vis_AxComput
172 (
173 DArray<Real> v_y,
174 DArray<Coord> velBuf,
175 DArray<Coord> position,
176 DArray<Real> alpha,
177 DArray<Attribute> attribute,
178 DArrayList<int> neighbor,
179 Real rest_density,
180 Real h,
181 Real dt,
182 DArray<Real> vis
183 )
184 {
185 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
186 if (pId >= position.size()) return;
187 if (!attribute[pId].isDynamic()) return;
188
189 Real tempValue = 0;
190 Coord Avi(0);
191 Real invAlpha_i = 1.0f / alpha[pId];
192
193 List<int>& list_i = neighbor[pId];
194 int nbSize = list_i.size();
195 for (int ne = 0; ne < nbSize; ne++)
196 {
197 int j = list_i[ne];
198 Real r = (position[pId] - position[j]).norm();
199 if (r > EPSILON)
200 {
201 Real wrr_ij = kernWRR(r, h);
202 Real ai = 1.0f / alpha[pId];
203 Real aj = 1.0f / alpha[j];
204 Coord nij = (position[j] - position[pId]) / r;
205 Coord vij = velBuf[pId] - velBuf[j];
206 tempValue = 2.5 * dt * (vis[pId] + vis[j]) * (ai + aj) * wrr_ij / rest_density;
207 Avi += tempValue * vij.dot(nij) * nij;
208 }
209 }
210 Avi += velBuf[pId];
211 v_y[3 * pId] = Avi[0];
212 v_y[3 * pId + 1] = Avi[1];
213 v_y[3 * pId + 2] = Avi[2];
214 }
215
216
217 template <typename Real, typename Coord>
218 __global__ void VC_Vis_r_Comput
219 (
220 DArray<Real> v_r,
221 DArray<Real> v_y,
222 DArray<Coord> vel_old,
223 DArray<Coord> position,
224 DArray<Attribute> attribute
225 )
226 {
227 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
228 if (pId >= position.size()) return;
229 if (!attribute[pId].isDynamic()) return;
230 Coord temp_vel = vel_old[pId];
231 v_r[3 * pId] = temp_vel[0] - v_y[3 * pId];
232 v_r[3 * pId + 1] = temp_vel[1] - v_y[3 * pId + 1];
233 v_r[3 * pId + 2] = temp_vel[2] - v_y[3 * pId + 2];
234 }
235
236 template <typename Real, typename Coord>
237 __global__ void VC_Vis_pToVector
238 (
239 DArray<Real> v_p,
240 DArray<Coord> pv,
241 DArray<Attribute> attribute
242 )
243 {
244 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
245 if (pId >= attribute.size()) return;
246 if (!attribute[pId].isDynamic()) return;
247 pv[pId][0] = v_p[3 * pId];
248 pv[pId][1] = v_p[3 * pId + 1];
249 pv[pId][2] = v_p[3 * pId + 2];
250 }
251
252
253 template <typename Real, typename Coord>
254 __global__ void VC_Vis_CoordToReal
255 (
256 DArray<Real> veloReal,
257 DArray<Coord> vel,
258 DArray<Attribute> attribute
259 )
260 {
261 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
262 if (pId >= attribute.size()) return;
263 if (!attribute[pId].isDynamic()) return;
264 veloReal[3 * pId] = vel[pId][0];
265 veloReal[3 * pId + 1] = vel[pId][1];
266 veloReal[3 * pId + 2] = vel[pId][2];
267 }
268
269 template <typename Real, typename Coord>
270 __global__ void VC_Vis_RealToVeloctiy
271 (
272 DArray<Coord> vel,
273 DArray<Real> veloReal,
274 DArray<Attribute> attribute
275 )
276 {
277 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
278 if (pId >= vel.size()) return;
279 if (!attribute[pId].isDynamic()) return;
280 vel[pId][0] = veloReal[3 * pId];
281 vel[pId][1] = veloReal[3 * pId + 1];
282 vel[pId][2] = veloReal[3 * pId + 2];
283 }
284
285 template <typename Real, typename Coord>
286 __global__ void VC_CrossVis
287 (
288 DArray<Real> crossVis,
289 DArray<Coord> velBuf,
290 DArray<Coord> position,
291 DArray<Real> alpha,
292 DArray<Attribute> attribute,
293 DArrayList<int> neighbor,
294 Real smoothingLength,
295 Real visTop,
296 Real visFloor,
297 Real K,
298 Real N
299 )
300 {
301 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
302 if (pId >= position.size()) return;
303 if (!attribute[pId].isDynamic()) return;
304
305 Coord pos_i = position[pId];
306 Coord vel_i = velBuf[pId];
307 Real invAlpha_i = 1.0f / alpha[pId];
308
309 Coord dv(0);
310 Coord vij(0);
311
312 List<int>& list_i = neighbor[pId];
313 int nbSize = list_i.size();
314 for (int ne = 0; ne < nbSize; ne++)
315 {
316 int j = list_i[ne];
317 Real r = (position[pId] - position[j]).norm();
318 if(r > EPSILON)
319 {
320 Real wr_ij = kernWR(r, smoothingLength);
321 Coord g = -invAlpha_i*(pos_i - position[j])*wr_ij*(1.0f / r);
322 vij = 0.5*(velBuf[pId] - velBuf[j]);
323 dv[0] += vij[0] * g[0];
324 dv[1] += vij[1] * g[1];
325 dv[2] += vij[2] * g[2];
326 }
327
328 }
329 Real Norm = VB_FNorm(dv[0], dv[1], dv[2]);
330 crossVis[pId] = VB_Viscosity(visTop, visFloor, Norm, K, N);
331 if (pId == 100) printf("viscosity : %f\r\n", crossVis[pId]);
332 }
333
334 template <typename Real, typename Coord>
335 __global__ void VC_updateScar
336 (
337 DArray<Real> scar,
338 DArray<Coord> velocity
339 )
340 {
341 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
342 if (pId >= velocity.size()) return;
343
344 scar[pId] = velocity[pId].norm();
345
346 }
347
348 __global__ void VC_AttributeInit
349 (
350 DArray<Attribute> atts
351 )
352 {
353 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
354 if (pId >= atts.size()) return;
355
356 atts[pId].setFluid();
357 atts[pId].setDynamic();
358 }
359
360 template <typename Real>
361 __global__ void VC_ViscosityValueUpdate
362 (
363 DArray<Real> viscosities,
364 Real viscosityValue
365 )
366 {
367 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
368 if (pId >= viscosities.size()) return;
369 viscosities[pId] = viscosityValue;
370 }
371
372 template<typename TDataType>
373 ApproximateImplicitViscosity<TDataType>::~ApproximateImplicitViscosity()
374 {
375 m_alpha.clear();
376 //if (m_reduce)
377 //{
378 // delete m_reduce;
379 //}
380 };
381
382
383 template<typename TDataType>
384 bool ApproximateImplicitViscosity<TDataType>::SetCross()
385 {
386 CrossVisCeil = this->varViscosity()->getValue();
387 CrossVisFloor = this->varLowerBoundViscosity()->getValue();
388 CrossVisFloor = CrossVisFloor < CrossVisCeil ? CrossVisFloor : CrossVisCeil;
389
390 Cross_K = this->varCrossK()->getValue();
391 Cross_N = this->varCrossN()->getValue();
392 std::cout << "*Non-Newtonian Fluid, Viscosity:" << CrossVisFloor
393 << " to " << CrossVisCeil << ", Cross_K:"
394 << Cross_K << ", Cross_N:" << Cross_N
395 << std::endl;
396 return true;
397 };
398
399 template<typename TDataType>
400 void ApproximateImplicitViscosity<TDataType>::constrain()
401 {
402 std::cout << "*Approximate Vicosity Solver::Dynamic Viscosity: " << this->varViscosity()->getValue() << std::endl;
403
404 int num = this->inPosition()->size();
405
406 if ((num != m_alpha.size())||
407 (num != velOld.size()) ||
408 (num != m_viscosity.size()))
409 {
410 m_alpha.resize(num);
411 v_y.resize(3 * num);
412 v_r.resize(3 * num);
413 v_p.resize(3 * num);
414 v_pv.resize(num);
415 m_VelocityReal.resize(3 * num);
416 velOld.resize(num);
417 velBuf.resize(num);
418 m_viscosity.resize(num);
419 m_reduce = Reduction<float>::Create(num);
420 m_arithmetic_v = Arithmetic<float>::Create(3 * num);
421 }
422
423 Real dt = this->inTimeStep()->getValue();
424
425 if (this->inAttribute()->isEmpty() || this->inAttribute()->size()!= this->inPosition()->size())
426 {
427 this->inAttribute()->allocate();
428 this->inAttribute()->resize(this->inPosition()->size());
429 cuExecute(num, VC_AttributeInit, this->inAttribute()->getData());
430 }
431
432 auto& m_position = this->inPosition()->getData();
433 auto& m_velocity = this->inVelocity()->getData();
434 auto& m_neighborhood = this->inNeighborIds()->getData();
435 auto& m_attribute = this->inAttribute()->getData();
436 auto m_smoothingLength = this->varSmoothingLength()->getValue();
437 auto m_restDensity = this->varRestDensity()->getValue();
438
439 cuExecute(num, VC_ViscosityValueUpdate,
440 m_viscosity,
441 this->varViscosity()->getValue()
442 );
443
444 m_alpha.reset();
445 cuExecute(num, VC_ComputeAlpha,
446 m_alpha,
447 m_position,
448 m_attribute,
449 m_neighborhood,
450 m_smoothingLength);
451
452 //VC_CorrectAlpha << <pDims, BLOCK_SIZE >> > (
453 // m_alpha,
454 // m_maxAlpha);
455
456 //IF Cross Model is active, viscous coefficients should be computed.
457
458
459 if (this->varFluidType()->getValue() == FluidType::NonNewtonianFluid)
460 {
461 SetCross();
462 cuExecute(num, VC_CrossVis,
463 m_viscosity,
464 m_velocity,
465 m_position,
466 m_alpha,
467 m_attribute,
468 m_neighborhood,
469 m_smoothingLength,
470 CrossVisCeil,
471 CrossVisFloor,
472 Cross_K,
473 Cross_N
474 );
475 }
476
477 velOld.assign(m_velocity);
478 v_y.reset();
479
480 cuExecute(num, VC_Vis_AxComput,
481 v_y,
482 m_velocity,
483 m_position,
484 m_alpha,
485 m_attribute,
486 m_neighborhood,
487 m_restDensity,
488 m_smoothingLength,
489 dt,
490 m_viscosity
491 );
492
493 v_r.reset();
494 cuExecute(num, VC_Vis_r_Comput,
495 v_r,
496 v_y,
497 velOld,
498 m_position,
499 m_attribute
500 );
501
502 v_p.assign(v_r);
503
504 Real Vrr = m_arithmetic_v->Dot(v_r, v_r);
505 Real Verr = sqrt(Vrr / v_r.size());
506 int VisItor = 0;
507 Real initErr = Verr;
508
509 std::cout <<"*Approximate Vicosity Solver::Residual:" << Vrr <<std::endl;
510 while (VisItor < 1000 && Verr / initErr > 0.01f && Vrr > 1000.0f * EPSILON)
511 {
512 VisItor++;
513 // //The type of "v_p" should convert to DArray<Coord>
514 cuExecute(num, VC_Vis_pToVector,
515 v_p,
516 v_pv,
517 m_attribute
518 );
519
520 v_y.reset();
521 cuExecute(num, VC_Vis_AxComput,
522 v_y,
523 v_pv,
524 m_position,
525 m_alpha,
526 m_attribute,
527 m_neighborhood,
528 m_restDensity,
529 m_smoothingLength,
530 dt,
531 m_viscosity
532 );
533
534 float alpha = Vrr / m_arithmetic_v->Dot(v_p, v_y);
535
536 cuExecute(num, VC_Vis_CoordToReal,
537 m_VelocityReal,
538 m_velocity,
539 m_attribute
540 );
541
542 Function2Pt::saxpy(m_VelocityReal, v_p, m_VelocityReal, alpha);
543
544 cuExecute(num, VC_Vis_RealToVeloctiy,
545 m_velocity,
546 m_VelocityReal,
547 m_attribute
548 );
549
550 Function2Pt::saxpy(v_r, v_y, v_r, -alpha);
551 Real Vrr_old = Vrr;
552
553 Vrr = m_arithmetic_v->Dot(v_r, v_r);
554 Real beta = Vrr / Vrr_old;
555 Function2Pt::saxpy(v_p, v_p, v_r, beta);
556
557 Verr = sqrt(Vrr / v_r.size());
558 }
559 std::cout << "*Approximate Vicosity Solver::Iteration #" << VisItor << ", Relative Resisual:" << Verr / initErr << std::endl;
560 };
561
562 DEFINE_CLASS(ApproximateImplicitViscosity);
563}
Generated by doxygen 1.13.2