Skip to content

Commit 484feb3

Browse files
neutrinocerosneutrinoceros
committed
ENH: implement Triangular Shaped Cloud deposition method
1 parent 89d7406 commit 484feb3

13 files changed

+305
-35
lines changed

‎gpgi/clib/_deposition_methods.pyx‎

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -81,3 +81,151 @@ def _deposit_pic_3D(
8181
j = hci_v[ipart][1]
8282
k = hci_v[ipart][2]
8383
out_v[i][j][k] += field_v[ipart]
84+
85+
86+
@cython.cdivision(True)
87+
@cython.boundscheck(False)
88+
@cython.wraparound(False)
89+
def _deposit_tsc_1D(
90+
np.ndarray[real, ndim=1] cell_edges_x1,
91+
np.ndarray[real, ndim=1] cell_edges_x2,
92+
np.ndarray[real, ndim=1] cell_edges_x3,
93+
np.ndarray[real, ndim=2] particle_coords,
94+
np.ndarray[real, ndim=1] field,
95+
np.ndarray[np.uint16_t, ndim=2] hci,
96+
np.ndarray[real, ndim=1] out,
97+
):
98+
cdef Py_ssize_t ipart, particle_count, i, oci
99+
particle_count = hci.shape[0]
100+
101+
cdef real x, d
102+
cdef np.uint16_t ci
103+
104+
cdef real[:] field_v = field
105+
cdef real[:] out_v = out
106+
107+
# weight arrays
108+
cdef real[3] w
109+
110+
for ipart in range(particle_count):
111+
x = particle_coords[ipart, 0]
112+
ci = hci[ipart, 0]
113+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
114+
w[0] = 0.5 * (1 - d) ** 2
115+
w[1] = 0.75 - (d - 0.5) ** 2
116+
w[2] = 0.5 * d**2
117+
118+
for i in range(3):
119+
oci = ci - 1 + i
120+
out_v[oci] += w[i] * field_v[ipart]
121+
122+
@cython.cdivision(True)
123+
@cython.boundscheck(False)
124+
@cython.wraparound(False)
125+
def _deposit_tsc_2D(
126+
np.ndarray[real, ndim=1] cell_edges_x1,
127+
np.ndarray[real, ndim=1] cell_edges_x2,
128+
np.ndarray[real, ndim=1] cell_edges_x3,
129+
np.ndarray[real, ndim=2] particle_coords,
130+
np.ndarray[real, ndim=1] field,
131+
np.ndarray[np.uint16_t, ndim=2] hci,
132+
np.ndarray[real, ndim=2] out,
133+
):
134+
cdef Py_ssize_t ipart, particle_count, i, j, oci, ocj
135+
particle_count = hci.shape[0]
136+
137+
cdef real x, d
138+
cdef np.uint16_t ci, cj
139+
140+
cdef real[:] field_v = field
141+
cdef real[:, :] out_v = out
142+
143+
# weight arrays
144+
cdef real[3] w1, w2
145+
cdef real[3][3] w
146+
147+
for ipart in range(particle_count):
148+
x = particle_coords[ipart, 0]
149+
ci = hci[ipart, 0]
150+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
151+
w1[0] = 0.5 * (1 - d) ** 2
152+
w1[1] = 0.75 - (d - 0.5) ** 2
153+
w1[2] = 0.5 * d**2
154+
155+
x = particle_coords[ipart, 1]
156+
cj = hci[ipart, 1]
157+
d = (x - cell_edges_x2[cj]) / (cell_edges_x2[cj + 1] - cell_edges_x2[cj])
158+
w2[0] = 0.5 * (1 - d) ** 2
159+
w2[1] = 0.75 - (d - 0.5) ** 2
160+
w2[2] = 0.5 * d**2
161+
162+
for i in range(3):
163+
for j in range(3):
164+
w[i][j] = w1[i] * w2[j]
165+
166+
for i in range(3):
167+
oci = ci - 1 + i
168+
for j in range(3):
169+
ocj = cj - 1 + j
170+
out_v[oci, ocj] += w[i][j] * field_v[ipart]
171+
172+
173+
@cython.cdivision(True)
174+
@cython.boundscheck(False)
175+
@cython.wraparound(False)
176+
def _deposit_tsc_3D(
177+
np.ndarray[real, ndim=1] cell_edges_x1,
178+
np.ndarray[real, ndim=1] cell_edges_x2,
179+
np.ndarray[real, ndim=1] cell_edges_x3,
180+
np.ndarray[real, ndim=2] particle_coords,
181+
np.ndarray[real, ndim=1] field,
182+
np.ndarray[np.uint16_t, ndim=2] hci,
183+
np.ndarray[real, ndim=3] out,
184+
):
185+
cdef Py_ssize_t ipart, particle_count, i, j, k, oci, ocj, ock
186+
particle_count = hci.shape[0]
187+
188+
cdef real x, d
189+
cdef np.uint16_t ci, cj, ck
190+
191+
cdef real[:] field_v = field
192+
cdef real[:, :, :] out_v = out
193+
194+
# weight arrays
195+
cdef real[3] w1, w2, w3
196+
cdef real[3][3][3] w
197+
198+
for ipart in range(particle_count):
199+
x = particle_coords[ipart, 0]
200+
ci = hci[ipart, 0]
201+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
202+
w1[0] = 0.5 * (1 - d) ** 2
203+
w1[1] = 0.75 - (d - 0.5) ** 2
204+
w1[2] = 0.5 * d**2
205+
206+
x = particle_coords[ipart, 1]
207+
cj = hci[ipart, 1]
208+
d = (x - cell_edges_x2[cj]) / (cell_edges_x2[cj + 1] - cell_edges_x2[cj])
209+
w2[0] = 0.5 * (1 - d) ** 2
210+
w2[1] = 0.75 - (d - 0.5) ** 2
211+
w2[2] = 0.5 * d**2
212+
213+
x = particle_coords[ipart, 2]
214+
ck = hci[ipart, 2]
215+
d = (x - cell_edges_x3[ck]) / (cell_edges_x3[ck + 1] - cell_edges_x3[ck])
216+
w3[0] = 0.5 * (1 - d) ** 2
217+
w3[1] = 0.75 - (d - 0.5) ** 2
218+
w3[2] = 0.5 * d**2
219+
220+
for i in range(3):
221+
for j in range(3):
222+
for k in range(3):
223+
w[i][j][k] = w1[i] * w2[j] * w3[k]
224+
225+
for i in range(3):
226+
oci = ci - 1 + i
227+
for j in range(3):
228+
ocj = cj - 1 + j
229+
for k in range(3):
230+
ock = ck - 1 + k
231+
out_v[oci, ocj, ock] += w[i][j][k] * field_v[ipart]

‎gpgi/lib/_deposition_methods.py‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,121 @@
1+
import numpy as np
2+
3+
14
def _deposit_pic(
25
cell_edges_x1, cell_edges_x2, cell_edges_x3, particle_coords, field, hci, out
36
):
47
for ipart in range(len(hci)):
58
md_idx = tuple(hci[ipart])
69
out[md_idx] += field[ipart]
10+
11+
12+
def _deposit_cic(
13+
cell_edges_x1, cell_edges_x2, cell_edges_x3, particle_coords, field, hci, out
14+
):
15+
return NotImplementedError("Cloud-in-cell deposition method is not implemented yet")
16+
17+
18+
def _deposit_tsc_1D(
19+
cell_edges_x1, cell_edges_x2, cell_edges_x3, particle_coords, field, hci, out
20+
):
21+
nparticles = hci.shape[0]
22+
23+
# weight array
24+
w = np.zeros(3, dtype=field.dtype)
25+
26+
for ipart in range(nparticles):
27+
x = particle_coords[ipart, 0]
28+
29+
ci = hci[ipart, 0]
30+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
31+
w[0] = 0.5 * (1 - d) ** 2
32+
w[1] = 0.75 - (d - 0.5) ** 2
33+
w[2] = 0.5 * d**2
34+
35+
for k, j in enumerate((ci - 1, ci, ci + 1)):
36+
out[j] += w[k] * field[ipart]
37+
38+
39+
def _deposit_tsc_2D(
40+
cell_edges_x1, cell_edges_x2, cell_edges_x3, particle_coords, field, hci, out
41+
):
42+
nparticles = hci.shape[0]
43+
44+
# weight arrays
45+
w1 = np.zeros(3, dtype=field.dtype)
46+
w2 = np.zeros(3, dtype=field.dtype)
47+
48+
w = np.zeros((3, 3), dtype=field.dtype)
49+
50+
for ipart in range(nparticles):
51+
x = particle_coords[ipart, 0]
52+
ci = hci[ipart, 0]
53+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
54+
assert d > 0
55+
w1[0] = 0.5 * (1 - d) ** 2
56+
w1[1] = 0.75 - (d - 0.5) ** 2
57+
w1[2] = 0.5 * d**2
58+
59+
x = particle_coords[ipart, 1]
60+
cj = hci[ipart, 1]
61+
d = (x - cell_edges_x2[cj]) / (cell_edges_x2[cj + 1] - cell_edges_x2[cj])
62+
assert d > 0
63+
w2[0] = 0.5 * (1 - d) ** 2
64+
w2[1] = 0.75 - (d - 0.5) ** 2
65+
w2[2] = 0.5 * d**2
66+
67+
for i in range(3):
68+
for j in range(3):
69+
w[i, j] = w1[i] * w2[j]
70+
71+
for i, oci in enumerate((ci - 1, ci, ci + 1)):
72+
for j, ocj in enumerate((cj - 1, cj, cj + 1)):
73+
out[oci, ocj] += w[i, j] * field[ipart]
74+
75+
76+
def _deposit_tsc_3D(
77+
cell_edges_x1, cell_edges_x2, cell_edges_x3, particle_coords, field, hci, out
78+
):
79+
nparticles = hci.shape[0]
80+
81+
# weight arrays
82+
w1 = np.zeros(3, dtype=field.dtype)
83+
w2 = np.zeros(3, dtype=field.dtype)
84+
w3 = np.zeros(3, dtype=field.dtype)
85+
86+
w = np.zeros((3, 3, 3), dtype=field.dtype)
87+
88+
for ipart in range(nparticles):
89+
x = particle_coords[ipart, 0]
90+
ci = hci[ipart, 0]
91+
d = (x - cell_edges_x1[ci]) / (cell_edges_x1[ci + 1] - cell_edges_x1[ci])
92+
assert d > 0
93+
w1[0] = 0.5 * (1 - d) ** 2
94+
w1[1] = 0.75 - (d - 0.5) ** 2
95+
w1[2] = 0.5 * d**2
96+
97+
x = particle_coords[ipart, 1]
98+
cj = hci[ipart, 1]
99+
d = (x - cell_edges_x2[cj]) / (cell_edges_x2[cj + 1] - cell_edges_x2[cj])
100+
assert d > 0
101+
w2[0] = 0.5 * (1 - d) ** 2
102+
w2[1] = 0.75 - (d - 0.5) ** 2
103+
w2[2] = 0.5 * d**2
104+
105+
x = particle_coords[ipart, 2]
106+
ck = hci[ipart, 2]
107+
d = (x - cell_edges_x3[ck]) / (cell_edges_x3[ck + 1] - cell_edges_x3[ck])
108+
assert d > 0
109+
w3[0] = 0.5 * (1 - d) ** 2
110+
w3[1] = 0.75 - (d - 0.5) ** 2
111+
w3[2] = 0.5 * d**2
112+
113+
for i in range(3):
114+
for j in range(3):
115+
for k in range(3):
116+
w[i, j, k] = w1[i] * w2[j] * w3[k]
117+
118+
for i, oci in enumerate((ci - 1, ci, ci + 1)):
119+
for j, ocj in enumerate((cj - 1, cj, cj + 1)):
120+
for k, ock in enumerate((ck - 1, ck, ck + 1)):
121+
out[oci, ocj, ock] += w[i, j, k] * field[ipart]

‎gpgi/types.py‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -333,6 +333,9 @@ def deposit(
333333
from .clib._deposition_methods import _deposit_pic_1D # type: ignore [import]
334334
from .clib._deposition_methods import _deposit_pic_2D # type: ignore [import]
335335
from .clib._deposition_methods import _deposit_pic_3D # type: ignore [import]
336+
from .clib._deposition_methods import _deposit_tsc_1D # type: ignore [import]
337+
from .clib._deposition_methods import _deposit_tsc_2D # type: ignore [import]
338+
from .clib._deposition_methods import _deposit_tsc_3D # type: ignore [import]
336339

337340
if not hasattr(self, "_cache"):
338341
self._cache: dict[tuple[Name, DepositionMethod], np.ndarray] = {}
@@ -364,7 +367,12 @@ def deposit(
364367
_deposit_pic_1D,
365368
_deposit_pic_2D,
366369
_deposit_pic_3D,
367-
]
370+
],
371+
DepositionMethod.TRIANGULAR_SHAPED_CLOUD: [
372+
_deposit_tsc_1D,
373+
_deposit_tsc_2D,
374+
_deposit_tsc_3D,
375+
],
368376
}
369377
if mkey not in known_methods:
370378
raise NotImplementedError(f"method {method} is not implemented yet")

‎tests/pytest_mpl_baseline/test_1D_deposit_geometric-pic.png‎

21 KB
Loading

‎tests/pytest_mpl_baseline/test_1D_deposit_geometric-tsc.png‎

21.3 KB
Loading

‎tests/pytest_mpl_baseline/test_1D_deposit_geometric.png‎

-17.5 KB
Binary file not shown.

‎tests/pytest_mpl_baseline/test_1D_deposit_linear-pic.png‎

20.8 KB
Loading

‎tests/pytest_mpl_baseline/test_1D_deposit_linear-tsc.png‎

21.3 KB
Loading

‎tests/pytest_mpl_baseline/test_1D_deposit_linear.png‎

-17.3 KB
Binary file not shown.

‎tests/pytest_mpl_baseline/test_2D_deposit_pic.png‎

33.6 KB
Loading

0 commit comments

Comments
 (0)