@@ -35,8 +35,8 @@ This function will setup the particle to node and particle to cell index for 2D
3535 Δdx = mp. ξ[ix, 1 ] - grid. ξ[mp. p2n[ix, iy], 1 ]
3636 Δdy = mp. ξ[ix, 2 ] - grid. ξ[mp. p2n[ix, iy], 2 ]
3737 # compute basis function
38- Nx, dNx = linearBasis (Δdx, grid. dx)
39- Ny, dNy = linearBasis (Δdy, grid. dy)
38+ Nx, dNx = linearbasis (Δdx, grid. dx)
39+ Ny, dNy = linearbasis (Δdy, grid. dy)
4040 mp. Nij[ix, iy] = Nx * Ny # shape function
4141 mp.∂Nx[ix, iy] = dNx * Ny # x-gradient shape function
4242 mp.∂Ny[ix, iy] = dNy * Nx # y-gradient shape function
@@ -70,9 +70,9 @@ This function will setup the particle to node and particle to cell index for 3D
7070 Δdy = mp. ξ[ix, 2 ] - grid. ξ[mp. p2n[ix, iy], 2 ]
7171 Δdz = mp. ξ[ix, 3 ] - grid. ξ[mp. p2n[ix, iy], 3 ]
7272 # compute basis function
73- Nx, dNx = linearBasis (Δdx, grid. dx)
74- Ny, dNy = linearBasis (Δdy, grid. dy)
75- Nz, dNz = linearBasis (Δdz, grid. dz)
73+ Nx, dNx = linearbasis (Δdx, grid. dx)
74+ Ny, dNy = linearbasis (Δdy, grid. dy)
75+ Nz, dNz = linearbasis (Δdz, grid. dz)
7676 mp. Nij[ix, iy] = Nx * Ny * Nz
7777 mp.∂Nx[ix, iy] = dNx * Ny * Nz # x-gradient shape function
7878 mp.∂Ny[ix, iy] = dNy * Nx * Nz # y-gradient shape function
@@ -181,7 +181,176 @@ end
181181@kernel inbounds = true function initmpstatus! (
182182 grid:: DeviceGrid2D{T1, T2} ,
183183 mp :: DeviceParticle2D{T1, T2} ,
184- :: Val{:bspline}
184+ :: Val{:bspline2}
185+ ) where {T1, T2}
186+ ix = @index (Global)
187+ if ix ≤ mp. np
188+ gdx_1, gdy_1 = inv (grid. dx), inv (grid. dy)
189+ # update mass and momentum
190+ mp. ms[ix] = mp. Ω[ix] * mp. ρs[ix]
191+ mp. ps[ix, 1 ] = mp. ms[ix] * mp. vs[ix, 1 ]
192+ mp. ps[ix, 2 ] = mp. ms[ix] * mp. vs[ix, 2 ]
193+ # base index in the grid
194+ # note the base index needs a shift of 0.5 (see Taichi)
195+ mξx, mξy = mp. ξ[ix, 1 ], mp. ξ[ix, 2 ]
196+ bnx = unsafe_trunc (T1, fld (mξx - T2 (0.5 ) * grid. dx - grid. x1, grid. dx))
197+ bny = unsafe_trunc (T1, fld (mξy - T2 (0.5 ) * grid. dy - grid. y1, grid. dy))
198+ bid = T1 (grid. nny * bnx + grid. nny - bny)
199+ gξx, gξy = grid. ξ[bid, 1 ], grid. ξ[bid, 2 ]
200+ # p2n index
201+ @KAunroll for iy in 1 : 9
202+ # p2n index
203+ bac = T1 (cld (iy, 3 ) - 1 )
204+ iyc = T1 (iy - bac * 3 )
205+ mp. p2n[ix, iy] = T1 (bid + bac * grid. nny - (iyc - 1 ))
206+ end
207+ # compute x value in the basis function N(x)
208+ x1 = gdx_1 * (mξx - gξx)
209+ x2 = gdx_1 * (mξx - gξx - grid. dx)
210+ x3 = gdx_1 * (mξx - gξx - grid. dx - grid. dx)
211+ y1 = gdy_1 * (mξy - gξy)
212+ y2 = gdy_1 * (mξy - gξy - grid. dy)
213+ y3 = gdy_1 * (mξy - gξy - grid. dy - grid. dy)
214+ Nx1, Nx2, Nx3, ∂x1, ∂x2, ∂x3 = bspline2basis (x1, x2, x3, gdx_1)
215+ Ny1, Ny2, Ny3, ∂y1, ∂y2, ∂y3 = bspline2basis (y1, y2, y3, gdy_1)
216+ # assign the value (in order)
217+ mp. Nij[ix, 1 ], mp. Nij[ix, 2 ], mp. Nij[ix, 3 ] = Nx1 * Ny1, Nx1 * Ny2, Nx1 * Ny3
218+ mp. Nij[ix, 4 ], mp. Nij[ix, 5 ], mp. Nij[ix, 6 ] = Nx2 * Ny1, Nx2 * Ny2, Nx2 * Ny3
219+ mp. Nij[ix, 7 ], mp. Nij[ix, 8 ], mp. Nij[ix, 9 ] = Nx3 * Ny1, Nx3 * Ny2, Nx3 * Ny3
220+ # assign the x-gradient (in order)
221+ mp.∂Nx[ix, 1 ], mp.∂Nx[ix, 2 ], mp.∂Nx[ix, 3 ] = ∂x1 * Ny1, ∂x1 * Ny2, ∂x1 * Ny3
222+ mp.∂Nx[ix, 4 ], mp.∂Nx[ix, 5 ], mp.∂Nx[ix, 6 ] = ∂x2 * Ny1, ∂x2 * Ny2, ∂x2 * Ny3
223+ mp.∂Nx[ix, 7 ], mp.∂Nx[ix, 8 ], mp.∂Nx[ix, 9 ] = ∂x3 * Ny1, ∂x3 * Ny2, ∂x3 * Ny3
224+ # assign the y-gradient (in order)
225+ mp.∂Ny[ix, 1 ], mp.∂Ny[ix, 2 ], mp.∂Ny[ix, 3 ] = Nx1 * ∂y1, Nx1 * ∂y2, Nx1 * ∂y3
226+ mp.∂Ny[ix, 4 ], mp.∂Ny[ix, 5 ], mp.∂Ny[ix, 6 ] = Nx2 * ∂y1, Nx2 * ∂y2, Nx2 * ∂y3
227+ mp.∂Ny[ix, 7 ], mp.∂Ny[ix, 8 ], mp.∂Ny[ix, 9 ] = Nx3 * ∂y1, Nx3 * ∂y2, Nx3 * ∂y3
228+ end
229+ end
230+
231+ @kernel inbounds = true function initmpstatus!! (
232+ grid:: DeviceGrid2D{T1, T2} ,
233+ mp :: DeviceParticle2D{T1, T2} ,
234+ :: Val{:bspline2}
235+ ) where {T1, T2}
236+ ix = @index (Global)
237+ if ix ≤ mp. np
238+ gdx_1, gdy_1 = inv (grid. dx), inv (grid. dy)
239+ # update mass and momentum
240+ mp. ms[ix] = mp. Ω[ix] * mp. ρs[ix]
241+ mp. ps[ix, 1 ] = mp. ms[ix] * mp. vs[ix, 1 ]
242+ mp. ps[ix, 2 ] = mp. ms[ix] * mp. vs[ix, 2 ]
243+ # base index in the grid
244+ # note the base index needs a shift of 0.5 (see Taichi)
245+ mξx, mξy = mp. ξ[ix, 1 ], mp. ξ[ix, 2 ]
246+ bnx = unsafe_trunc (T1, fld (mξx - T2 (0.5 ) * grid. dx - grid. x1, grid. dx))
247+ bny = unsafe_trunc (T1, fld (mξy - T2 (0.5 ) * grid. dy - grid. y1, grid. dy))
248+ bid = T1 (grid. nny * bnx + grid. nny - bny)
249+ gξx, gξy = grid. ξ[bid, 1 ], grid. ξ[bid, 2 ]
250+ # p2n index
251+ @KAunroll for iy in 1 : 9
252+ # p2n index
253+ bac = T1 (cld (iy, 3 ) - 1 )
254+ iyc = T1 (iy - bac * 3 )
255+ mp. p2n[ix, iy] = T1 (bid + bac * grid. nny - (iyc - 1 ))
256+ end
257+ # compute x value in the basis function N(x)
258+ x1 = gdx_1 * (mξx - gξx)
259+ x2 = gdx_1 * (mξx - gξx - grid. dx)
260+ x3 = gdx_1 * (mξx - gξx - grid. dx - grid. dx)
261+ y1 = gdy_1 * (mξy - gξy)
262+ y2 = gdy_1 * (mξy - gξy - grid. dy)
263+ y3 = gdy_1 * (mξy - gξy - grid. dy - grid. dy)
264+ Nx1, Nx2, Nx3, ∂x1, ∂x2, ∂x3 = bspline2basis (x1, x2, x3, gdx_1)
265+ Ny1, Ny2, Ny3, ∂y1, ∂y2, ∂y3 = bspline2basis (y1, y2, y3, gdy_1)
266+ # assign the value (in order)
267+ mp. Nij[ix, 1 ], mp. Nij[ix, 2 ], mp. Nij[ix, 3 ] = Nx1 * Ny1, Nx1 * Ny2, Nx1 * Ny3
268+ mp. Nij[ix, 4 ], mp. Nij[ix, 5 ], mp. Nij[ix, 6 ] = Nx2 * Ny1, Nx2 * Ny2, Nx2 * Ny3
269+ mp. Nij[ix, 7 ], mp. Nij[ix, 8 ], mp. Nij[ix, 9 ] = Nx3 * Ny1, Nx3 * Ny2, Nx3 * Ny3
270+ if size (mp.∂Nx, 1 ) == mp. np
271+ # assign the x-gradient (in order)
272+ mp.∂Nx[ix, 1 ], mp.∂Nx[ix, 2 ], mp.∂Nx[ix, 3 ] = ∂x1 * Ny1, ∂x1 * Ny2, ∂x1 * Ny3
273+ mp.∂Nx[ix, 4 ], mp.∂Nx[ix, 5 ], mp.∂Nx[ix, 6 ] = ∂x2 * Ny1, ∂x2 * Ny2, ∂x2 * Ny3
274+ mp.∂Nx[ix, 7 ], mp.∂Nx[ix, 8 ], mp.∂Nx[ix, 9 ] = ∂x3 * Ny1, ∂x3 * Ny2, ∂x3 * Ny3
275+ # assign the y-gradient (in order)
276+ mp.∂Ny[ix, 1 ], mp.∂Ny[ix, 2 ], mp.∂Ny[ix, 3 ] = Nx1 * ∂y1, Nx1 * ∂y2, Nx1 * ∂y3
277+ mp.∂Ny[ix, 4 ], mp.∂Ny[ix, 5 ], mp.∂Ny[ix, 6 ] = Nx2 * ∂y1, Nx2 * ∂y2, Nx2 * ∂y3
278+ mp.∂Ny[ix, 7 ], mp.∂Ny[ix, 8 ], mp.∂Ny[ix, 9 ] = Nx3 * ∂y1, Nx3 * ∂y2, Nx3 * ∂y3
279+ end
280+ end
281+ end
282+
283+ @kernel inbounds = true function initmpstatus! (
284+ grid:: DeviceGrid3D{T1, T2} ,
285+ mp :: DeviceParticle3D{T1, T2} ,
286+ :: Val{:bspline2}
287+ ) where {T1, T2}
288+ ix = @index (Global)
289+ if ix ≤ mp. np
290+ gdx_1, gdy_1, gdz_1 = inv (grid. dx), inv (grid. dy), inv (grid. dz)
291+ # update mass and momentum
292+ mp. ms[ix] = mp. Ω[ix] * mp. ρs[ix]
293+ mp. ps[ix, 1 ] = mp. ms[ix] * mp. vs[ix, 1 ]
294+ mp. ps[ix, 2 ] = mp. ms[ix] * mp. vs[ix, 2 ]
295+ mp. ps[ix, 3 ] = mp. ms[ix] * mp. vs[ix, 3 ]
296+ # base index in the grid
297+ # note the base index needs a shift of 0.5 (see Taichi)
298+ mξx, mξy, mξz = mp. ξ[ix, 1 ], mp. ξ[ix, 2 ], mp. ξ[ix, 3 ]
299+ bnx = unsafe_trunc (T1, fld (mξx - T2 (0.5 ) * grid. dx - grid. x1, grid. dx))
300+ bny = unsafe_trunc (T1, fld (mξy - T2 (0.5 ) * grid. dy - grid. y1, grid. dy))
301+ bnz = unsafe_trunc (T1, fld (mξz - T2 (0.5 ) * grid. dz - grid. z1, grid. dz))
302+ bid = T1 (grid. nnx * grid. nny * bnz + grid. nny * bnx + bny + 1 )
303+ gξx, gξy, gξz = grid. ξ[bid, 1 ], grid. ξ[bid, 2 ], grid. ξ[bid, 3 ]
304+ # p2n index
305+ @KAunroll for k in 0 : 2
306+ for j in 0 : 2
307+ for i in 0 : 2
308+ iy = 1 + i + 3 * j + 9 * k
309+ mp. p2n[ix, iy] = bid + i + grid. nny * j + grid. nny * grid. nnx * k
310+ end
311+ end
312+ end
313+ # compute x value in the basis function N(x)
314+ x1 = gdx_1 * (mξx - gξx)
315+ x2 = gdx_1 * (mξx - gξx - grid. dx)
316+ x3 = gdx_1 * (mξx - gξx - grid. dx - grid. dx)
317+ y1 = gdy_1 * (mξy - gξy)
318+ y2 = gdy_1 * (mξy - gξy - grid. dy)
319+ y3 = gdy_1 * (mξy - gξy - grid. dy - grid. dy)
320+ z1 = gdz_1 * (mξz - gξz)
321+ z2 = gdz_1 * (mξz - gξz - grid. dz)
322+ z3 = gdz_1 * (mξz - gξz - grid. dz - grid. dz)
323+ Nx1, Nx2, Nx3, dx1, dx2, dx3 = bspline2basis (x1, x2, x3, gdx_1)
324+ Ny1, Ny2, Ny3, dy1, dy2, dy3 = bspline2basis (y1, y2, y3, gdy_1)
325+ Nz1, Nz2, Nz3, dz1, dz2, dz3 = bspline2basis (z1, z2, z3, gdz_1)
326+ # assign the value (in order)
327+ it = Int32 (1 )
328+ Nxs = (Nx1, Nx2, Nx3)
329+ Nys = (Ny1, Ny2, Ny3)
330+ Nzs = (Nz1, Nz2, Nz3)
331+ dxs = (dx1, dx2, dx3)
332+ dys = (dy1, dy2, dy3)
333+ dzs = (dz1, dz2, dz3)
334+ @KAunroll for k in 1 : 3 # z-direction
335+ for i in 1 : 3 # x-direction
336+ for j in 1 : 3 # y-direction
337+ mp. Nij[ix, it] = Nxs[i] * Nys[j] * Nzs[k]
338+ if size (mp.∂Nx, 1 ) == mp. np
339+ mp.∂Nx[ix, it] = dxs[i] * Nys[j] * Nzs[k]
340+ mp.∂Ny[ix, it] = Nxs[i] * dys[j] * Nzs[k]
341+ mp.∂Nz[ix, it] = Nxs[i] * Nys[j] * dzs[k]
342+ end
343+ it += Int32 (1 )
344+ end
345+ end
346+ end
347+ end
348+ end
349+
350+ @kernel inbounds = true function initmpstatus! (
351+ grid:: DeviceGrid2D{T1, T2} ,
352+ mp :: DeviceParticle2D{T1, T2} ,
353+ :: Val{:bspline3}
185354) where {T1, T2}
186355 ix = @index (Global)
187356 if ix ≤ mp. np
201370 ry = (mp. ξ[ix, 2 ] - grid. ξ[p2n, 2 ]) / grid. dy
202371 # compute basis function
203372 type_vx = get_type (grid. ξ[p2n, 1 ], grid. x1, grid. x2, grid. dx)
204- Nx, dNx = bsplinebasis (rx, grid. dx, type_vx)
373+ Nx, dNx = bspline3basis (rx, grid. dx, type_vx)
205374 type_vy = get_type (grid. ξ[p2n, 2 ], grid. y1, grid. y2, grid. dy)
206- Ny, dNy = bsplinebasis (ry, grid. dy, type_vy)
375+ Ny, dNy = bspline3basis (ry, grid. dy, type_vy)
207376 mp. Nij[ix, iy] = Nx * Ny
208377 mp.∂Nx[ix, iy] = dNx * Ny # x-gradient shape function
209378 mp.∂Ny[ix, iy] = dNy * Nx # y-gradient shape function
215384@kernel inbounds = true function initmpstatus! (
216385 grid:: DeviceGrid3D{T1, T2} ,
217386 mp :: DeviceParticle3D{T1, T2} ,
218- :: Val{:bspline }
387+ :: Val{:bspline3 }
219388) where {T1, T2}
220389 ix = @index (Global)
221390 if ix ≤ mp. np
@@ -242,11 +411,11 @@ end
242411 rz = (mpξ3 - grid. ξ[p2n, 3 ]) / grid. dz
243412 # compute basis function
244413 type_vx = get_type (grid. ξ[p2n, 1 ], grid. x1, grid. x2, grid. dx)
245- Nx, dNx = bsplinebasis (rx, grid. dx, type_vx)
414+ Nx, dNx = bspline3basis (rx, grid. dx, type_vx)
246415 type_vy = get_type (grid. ξ[p2n, 2 ], grid. y1, grid. y2, grid. dy)
247- Ny, dNy = bsplinebasis (ry, grid. dy, type_vy)
416+ Ny, dNy = bspline3basis (ry, grid. dy, type_vy)
248417 type_vz = get_type (grid. ξ[p2n, 3 ], grid. z1, grid. z2, grid. dz)
249- Nz, dNz = bsplinebasis (rz, grid. dz, type_vz)
418+ Nz, dNz = bspline3basis (rz, grid. dz, type_vz)
250419 mp. Nij[ix, iy] = T2 ( Nx * Ny * Nz)
251420 mp.∂Nx[ix, iy] = T2 (dNx * Ny * Nz) # x-gradient basis function
252421 mp.∂Ny[ix, iy] = T2 (dNy * Nx * Nz) # y-gradient basis function
0 commit comments