PeriDyno 1.0.0
Loading...
Searching...
No Matches
SemiImplicitDensitySolver.cu
Go to the documentation of this file.
1#include "SemiImplicitDensitySolver.h"
2
3#include "SummationDensity.h"
4
5namespace dyno
6{
7 template<typename TDataType>
8 SemiImplicitDensitySolver<TDataType>::SemiImplicitDensitySolver()
9 : ParticleApproximation<TDataType>()
10 {
11 this->varIterationNumber()->setValue(3);
12 this->varRestDensity()->setValue(Real(1000));
13
14 this->inAttribute()->tagOptional(true);
15
16 mSummation = std::make_shared<SummationDensity<TDataType>>();
17
18 this->inSmoothingLength()->connect(mSummation->inSmoothingLength());
19 this->inSamplingDistance()->connect(mSummation->inSamplingDistance());
20 this->inPosition()->connect(mSummation->inPosition());
21 this->inNeighborIds()->connect(mSummation->inNeighborIds());
22
23 mSummation->outDensity()->connect(this->outDensity());
24 }
25
26 template<typename TDataType>
27 SemiImplicitDensitySolver<TDataType>::~SemiImplicitDensitySolver()
28 {
29 mPosBuf.clear();
30 mPosOld.clear();
31 }
32
33#define CONSERVE_MOMETNUM
34#define CASE_0
35//#define CASE_1
36//#define CASE_2
37
38#define CHEBYSHEV_ACCELERATION
39//#define ANDERSON_ACCELERATION
40
41#ifdef CHEBYSHEV_ACCELERATION
42 template<typename Real, typename Coord>
43 __global__ void SSDS_Blend(
44 DArray<Coord> pos,
45 DArray<Coord> pos_k,
46 DArray<Coord> pos_k_minus,
47 Real omega)
48 {
49 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
50 if (pId >= pos.size()) return;
51
52 pos[pId] = omega * (pos_k[pId] - pos_k_minus[pId]) + pos_k_minus[pId];
53 }
54#endif
55
56#ifdef ANDERSON_ACCELERATION
57 template<typename Coord>
58 __global__ void SSDS_CalculateG(
59 DArray<Coord> g,
60 DArray<Coord> x,
61 DArray<Coord> fx)
62 {
63 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
64 if (pId >= g.size()) return;
65
66 g[pId] = fx[pId] - x[pId];
67 }
68
69 template<typename Real, typename Coord>
70 __global__ void SSDS_CalculateLSM(
71 DArray<Real> denom,
72 DArray<Real> numer,
73 DArray<Coord> g_k_minus,
74 DArray<Coord> g_k)
75 {
76 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
77 if (pId >= denom.size()) return;
78
79 Coord dg = g_k[pId] - g_k_minus[pId];
80
81 denom[pId] = dg.normSquared();
82 numer[pId] = dg.dot(g_k[pId]);
83 }
84
85 template<typename Real, typename Coord>
86 __global__ void SSDS_Blend(
87 DArray<Coord> pos,
88 DArray<Coord> pos_k,
89 DArray<Coord> pos_k_minus,
90 Real omega)
91 {
92 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
93 if (pId >= pos.size()) return;
94
95 pos[pId] = pos_k[pId] + omega * (pos_k_minus[pId] - pos_k[pId]);
96 }
97#endif
98
99
100 template <typename Real, typename Coord, typename Kernel>
101 __global__ void SSDS_OneJacobiStep(
102 DArray<Coord> posNext,
103 DArray<Coord> posBuf,
104 DArray<Coord> posOld,
105 DArray<Real> diagnals,
106 DArray<Real> rho,
107 DArrayList<int> neighbors,
108 Real smoothingLength,
109 Real samplingDistance,
110 Real kappa,
111 Real dt,
112 Kernel gradient,
113 Real scale)
114 {
115 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
116 if (pId >= posNext.size()) return;
117
118 Real rho_0 = Real(1000);
119
120 Real rho_i = rho[pId];
121 rho_i = rho_i > rho_0 ? rho_i : rho_0;
122
123 Real A = kappa * dt * dt / rho_0;
124 Real C_plus = rho_i / rho_0;
125 Real C_minus = Real(-1);
126
127 Coord pos_i = posBuf[pId];
128
129 List<int>& list_i = neighbors[pId];
130 int nbSize = list_i.size();
131
132 Real a_i = Real(1);
133 Coord posAcc_i = posOld[pId];
134 Coord deltaPos_i = Coord(0);
135 for (int ne = 0; ne < nbSize; ne++)
136 {
137 int j = list_i[ne];
138 Coord pos_j = posBuf[j];
139 Real r = (pos_i - pos_j).norm();
140
141 if (r > EPSILON)
142 {
143 Real a_ij = A * gradient(r, smoothingLength, scale) * (1.0f / r);
144
145#ifdef CONSERVE_MOMETNUM
146 posAcc_i += C_minus * a_ij * pos_j + C_plus * a_ij * (pos_j - pos_i);
147
148 Coord posAcc_ji = C_minus * a_ij * (pos_i) + C_plus * a_ij * (pos_i - pos_j);
149 Real a_ji = C_minus * a_ij;
150
151 atomicAdd(&posNext[j][0], posAcc_ji[0]);
152 atomicAdd(&posNext[j][1], posAcc_ji[1]);
153 atomicAdd(&posNext[j][2], posAcc_ji[2]);
154 atomicAdd(&diagnals[j], a_ji);
155#else
156 posAcc_i += C_minus * a_ij * pos_j + C_plus * a_ij * (pos_j - pos_i);
157#endif // CONSERVE_MOMETNUM
158 a_i += C_minus * a_ij;
159 }
160 }
161
162#ifdef CONSERVE_MOMETNUM
163 //To ensure momentum conservation
164 atomicAdd(&posNext[pId][0], posAcc_i[0]);
165 atomicAdd(&posNext[pId][1], posAcc_i[1]);
166 atomicAdd(&posNext[pId][2], posAcc_i[2]);
167
168 atomicAdd(&diagnals[pId], a_i);
169#else
170 posNext[pId] = posAcc_i / a_i;;
171#endif // CONSERVE_MOMETNUM
172 }
173
174 template <typename Real, typename Coord>
175 __global__ void SSDS_ComputeNewPos(
176 DArray<Coord> posNext,
177 DArray<Real> diagnals)
178 {
179 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
180 if (pId >= posNext.size()) return;
181
182 posNext[pId] = posNext[pId] / diagnals[pId];
183 }
184
185 template <typename Real, typename Coord>
186 __global__ void SSDS_ComputeNewPos(
187 DArray<Coord> posNext,
188 DArray<Coord> posOld,
189 DArray<Attribute> attributes,
190 DArray<Real> diagnals)
191 {
192 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
193 if (pId >= posNext.size()) return;
194
195 posNext[pId] = attributes[pId].isDynamic() ? posNext[pId] / diagnals[pId] : posOld[pId];
196 }
197
198 template<typename TDataType>
199 void SemiImplicitDensitySolver<TDataType>::compute()
200 {
201 updatePosition();
202
203 updateVelocity();
204 }
205
206
207 template<typename TDataType>
208 void SemiImplicitDensitySolver<TDataType>::updatePosition()
209 {
210 int num = this->inPosition()->size();
211 Real dt = this->inTimeStep()->getValue();
212 auto& inPos = this->inPosition()->getData();
213
214 if (mPosOld.size() != num)
215 mPosOld.resize(num);
216
217 if (mPosBuf.size() != num)
218 mPosBuf.resize(num);
219
220 if (mDiagnals.size() != num)
221 mDiagnals.resize(num);
222
223 mPosOld.assign(inPos);
224
225 int itNum = this->varIterationNumber()->getValue();
226 int it = 0;
227 while (it < itNum)
228 {
229 mSummation->varRestDensity()->setValue(this->varRestDensity()->getValue());
230 mSummation->varKernelType()->setCurrentKey(this->varKernelType()->currentKey());
231 mSummation->update();
232
233 mPosBuf.assign(inPos);
234#ifdef CONSERVE_MOMETNUM
235 inPos.reset();
236#endif
237 mDiagnals.reset();
238 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
239 SSDS_OneJacobiStep,
240 inPos,
241 mPosBuf,
242 mPosOld,
243 mDiagnals,
244 mSummation->outDensity()->getData(),
245 this->inNeighborIds()->getData(),
246 this->inSmoothingLength()->getValue(),
247 this->inSamplingDistance()->getValue(),
248 this->varKappa()->getValue(),
249 dt);
250
251#ifdef CONSERVE_MOMETNUM
252 if (this->inAttribute()->isEmpty())
253 {
254 cuExecute(num,
255 SSDS_ComputeNewPos,
256 inPos,
257 mDiagnals);
258 }
259 else
260 {
261 cuExecute(num,
262 SSDS_ComputeNewPos,
263 inPos,
264 mPosOld,
265 this->inAttribute()->constData(),
266 mDiagnals);
267 }
268
269#endif
270
271 it++;
272 }
273 }
274
275 template <typename Real, typename Coord>
276 __global__ void SSDS_UpdateVelocity(
277 DArray<Coord> velArr,
278 DArray<Coord> curPos,
279 DArray<Coord> prePos,
280 Real dt)
281 {
282 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
283 if (pId >= velArr.size()) return;
284
285 velArr[pId] += (curPos[pId] - prePos[pId]) / dt;
286 }
287
288 template<typename TDataType>
289 void SemiImplicitDensitySolver<TDataType>::updateVelocity()
290 {
291 int num = this->inPosition()->size();
292
293 Real dt = this->inTimeStep()->getData();
294
295 cuExecute(num, SSDS_UpdateVelocity,
296 this->inVelocity()->getData(),
297 this->inPosition()->getData(),
298 mPosOld,
299 dt);
300 }
301
302 DEFINE_CLASS(SemiImplicitDensitySolver);
303}
Generated by doxygen 1.13.2