1#include "SemiImplicitDensitySolver.h"
3#include "SummationDensity.h"
7 template<typename TDataType>
8 SemiImplicitDensitySolver<TDataType>::SemiImplicitDensitySolver()
9 : ParticleApproximation<TDataType>()
11 this->varIterationNumber()->setValue(3);
12 this->varRestDensity()->setValue(Real(1000));
14 this->inAttribute()->tagOptional(true);
16 mSummation = std::make_shared<SummationDensity<TDataType>>();
18 this->inSmoothingLength()->connect(mSummation->inSmoothingLength());
19 this->inSamplingDistance()->connect(mSummation->inSamplingDistance());
20 this->inPosition()->connect(mSummation->inPosition());
21 this->inNeighborIds()->connect(mSummation->inNeighborIds());
23 mSummation->outDensity()->connect(this->outDensity());
26 template<typename TDataType>
27 SemiImplicitDensitySolver<TDataType>::~SemiImplicitDensitySolver()
33#define CONSERVE_MOMETNUM
38#define CHEBYSHEV_ACCELERATION
39//#define ANDERSON_ACCELERATION
41#ifdef CHEBYSHEV_ACCELERATION
42 template<typename Real, typename Coord>
43 __global__ void SSDS_Blend(
46 DArray<Coord> pos_k_minus,
49 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
50 if (pId >= pos.size()) return;
52 pos[pId] = omega * (pos_k[pId] - pos_k_minus[pId]) + pos_k_minus[pId];
56#ifdef ANDERSON_ACCELERATION
57 template<typename Coord>
58 __global__ void SSDS_CalculateG(
63 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
64 if (pId >= g.size()) return;
66 g[pId] = fx[pId] - x[pId];
69 template<typename Real, typename Coord>
70 __global__ void SSDS_CalculateLSM(
73 DArray<Coord> g_k_minus,
76 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
77 if (pId >= denom.size()) return;
79 Coord dg = g_k[pId] - g_k_minus[pId];
81 denom[pId] = dg.normSquared();
82 numer[pId] = dg.dot(g_k[pId]);
85 template<typename Real, typename Coord>
86 __global__ void SSDS_Blend(
89 DArray<Coord> pos_k_minus,
92 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
93 if (pId >= pos.size()) return;
95 pos[pId] = pos_k[pId] + omega * (pos_k_minus[pId] - pos_k[pId]);
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,
107 DArrayList<int> neighbors,
108 Real smoothingLength,
109 Real samplingDistance,
115 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
116 if (pId >= posNext.size()) return;
118 Real rho_0 = Real(1000);
120 Real rho_i = rho[pId];
121 rho_i = rho_i > rho_0 ? rho_i : rho_0;
123 Real A = kappa * dt * dt / rho_0;
124 Real C_plus = rho_i / rho_0;
125 Real C_minus = Real(-1);
127 Coord pos_i = posBuf[pId];
129 List<int>& list_i = neighbors[pId];
130 int nbSize = list_i.size();
133 Coord posAcc_i = posOld[pId];
134 Coord deltaPos_i = Coord(0);
135 for (int ne = 0; ne < nbSize; ne++)
138 Coord pos_j = posBuf[j];
139 Real r = (pos_i - pos_j).norm();
143 Real a_ij = A * gradient(r, smoothingLength, scale) * (1.0f / r);
145#ifdef CONSERVE_MOMETNUM
146 posAcc_i += C_minus * a_ij * pos_j + C_plus * a_ij * (pos_j - pos_i);
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;
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);
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;
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]);
168 atomicAdd(&diagnals[pId], a_i);
170 posNext[pId] = posAcc_i / a_i;;
171#endif // CONSERVE_MOMETNUM
174 template <typename Real, typename Coord>
175 __global__ void SSDS_ComputeNewPos(
176 DArray<Coord> posNext,
177 DArray<Real> diagnals)
179 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
180 if (pId >= posNext.size()) return;
182 posNext[pId] = posNext[pId] / diagnals[pId];
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)
192 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
193 if (pId >= posNext.size()) return;
195 posNext[pId] = attributes[pId].isDynamic() ? posNext[pId] / diagnals[pId] : posOld[pId];
198 template<typename TDataType>
199 void SemiImplicitDensitySolver<TDataType>::compute()
207 template<typename TDataType>
208 void SemiImplicitDensitySolver<TDataType>::updatePosition()
210 int num = this->inPosition()->size();
211 Real dt = this->inTimeStep()->getValue();
212 auto& inPos = this->inPosition()->getData();
214 if (mPosOld.size() != num)
217 if (mPosBuf.size() != num)
220 if (mDiagnals.size() != num)
221 mDiagnals.resize(num);
223 mPosOld.assign(inPos);
225 int itNum = this->varIterationNumber()->getValue();
229 mSummation->varRestDensity()->setValue(this->varRestDensity()->getValue());
230 mSummation->varKernelType()->setCurrentKey(this->varKernelType()->currentKey());
231 mSummation->update();
233 mPosBuf.assign(inPos);
234#ifdef CONSERVE_MOMETNUM
238 cuFirstOrder(num, this->varKernelType()->getDataPtr()->currentKey(), this->mScalingFactor,
244 mSummation->outDensity()->getData(),
245 this->inNeighborIds()->getData(),
246 this->inSmoothingLength()->getValue(),
247 this->inSamplingDistance()->getValue(),
248 this->varKappa()->getValue(),
251#ifdef CONSERVE_MOMETNUM
252 if (this->inAttribute()->isEmpty())
265 this->inAttribute()->constData(),
275 template <typename Real, typename Coord>
276 __global__ void SSDS_UpdateVelocity(
277 DArray<Coord> velArr,
278 DArray<Coord> curPos,
279 DArray<Coord> prePos,
282 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
283 if (pId >= velArr.size()) return;
285 velArr[pId] += (curPos[pId] - prePos[pId]) / dt;
288 template<typename TDataType>
289 void SemiImplicitDensitySolver<TDataType>::updateVelocity()
291 int num = this->inPosition()->size();
293 Real dt = this->inTimeStep()->getData();
295 cuExecute(num, SSDS_UpdateVelocity,
296 this->inVelocity()->getData(),
297 this->inPosition()->getData(),
302 DEFINE_CLASS(SemiImplicitDensitySolver);