Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
GreenFuncGravCartesian.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
18#include "shambase/aliases_int.hpp"
19#include "shammath/symtensor_collections.hpp"
20
21namespace shamphys {
22
36 template<class T, u32 low_order, u32 high_order>
37 class GreenFuncGravCartesian {
38 public:
39 inline static shammath::SymTensorCollection<T, low_order, high_order> get_der_tensors(
40 const sycl::vec<T, 3> &r);
41 };
42
43 template<class T, u32 low_order, u32 high_order>
44 inline shammath::SymTensorCollection<T, low_order, high_order> green_func_grav_cartesian(
45 const sycl::vec<T, 3> &r) {
46 return GreenFuncGravCartesian<T, low_order, high_order>::get_der_tensors(r);
47 }
48
50// Implementations for all cases
51// -----------
52// Do not look if you donw want your eyes to bleed, this is very VERY ugly
54#ifndef DOXYGEN
55 template<class T>
56 class GreenFuncGravCartesian<T, 0, 5> {
57 public:
58 inline static shammath::SymTensorCollection<T, 0, 5> get_der_tensors(
59 const sycl::vec<T, 3> &r) {
60 using namespace shammath;
61
62 T r1 = r.x();
63 T r2 = r.y();
64 T r3 = r.z();
65
66 T r1pow2 = r.x() * r.x();
67 T r2pow2 = r.y() * r.y();
68 T r3pow2 = r.z() * r.z();
69
70 T r1pow3 = r.x() * r1pow2;
71 T r2pow3 = r.y() * r2pow2;
72 T r3pow3 = r.z() * r3pow2;
73
74 T r1pow4 = r.x() * r1pow3;
75 T r2pow4 = r.y() * r2pow3;
76 T r3pow4 = r.z() * r3pow3;
77
78 T r1pow5 = r.x() * r1pow4;
79 T r2pow5 = r.y() * r2pow4;
80 T r3pow5 = r.z() * r3pow4;
81
82 T rsq = r1pow2 + r2pow2 + r3pow2;
83
84 T rnorm = sycl::sqrt(rsq);
85
86 T rm2 = 1 / (rsq);
87
88 T g0 = 1 / rnorm;
89 T g1 = -1 * rm2 * g0;
90 T g2 = -3 * rm2 * g1;
91 T g3 = -5 * rm2 * g2;
92 T g4 = -7 * rm2 * g3;
93 T g5 = -9 * rm2 * g4;
94
95 auto D5 = SymTensor3d_5<T>{
96 15 * g3 * r1 + 10 * g4 * r1pow3 + g5 * r1pow5,
97 (3 * g3 + 6 * g4 * r1pow2 + g5 * r1pow4) * r2,
98 (3 * g3 + 6 * g4 * r1pow2 + g5 * r1pow4) * r3,
99 3 * g3 * r1 + g4 * r1pow3 + 3 * g4 * r1 * r2pow2 + g5 * r1pow3 * r2pow2,
100 (3 * g4 * r1 + g5 * r1pow3) * r2 * r3,
101 3 * g3 * r1 + g4 * r1pow3 + 3 * g4 * r1 * r3pow2 + g5 * r1pow3 * r3pow2,
102 3 * g3 * r2 + 3 * g4 * r1pow2 * r2 + g4 * r2pow3 + g5 * r1pow2 * r2pow3,
103 (g3 + g4 * r1pow2 + g4 * r2pow2 + g5 * r1pow2 * r2pow2) * r3,
104 r2 * (g3 + g4 * r1pow2 + g4 * r3pow2 + g5 * r1pow2 * r3pow2),
105 3 * g3 * r3 + 3 * g4 * r1pow2 * r3 + g4 * r3pow3 + g5 * r1pow2 * r3pow3,
106 r1 * (3 * g3 + 6 * g4 * r2pow2 + g5 * r2pow4),
107 r1 * (3 * g4 * r2 + g5 * r2pow3) * r3,
108 r1 * (g3 + g4 * r2pow2 + g4 * r3pow2 + g5 * r2pow2 * r3pow2),
109 r1 * r2 * (3 * g4 * r3 + g5 * r3pow3),
110 r1 * (3 * g3 + 6 * g4 * r3pow2 + g5 * r3pow4),
111 15 * g3 * r2 + 10 * g4 * r2pow3 + g5 * r2pow5,
112 (3 * g3 + 6 * g4 * r2pow2 + g5 * r2pow4) * r3,
113 3 * g3 * r2 + g4 * r2pow3 + 3 * g4 * r2 * r3pow2 + g5 * r2pow3 * r3pow2,
114 3 * g3 * r3 + 3 * g4 * r2pow2 * r3 + g4 * r3pow3 + g5 * r2pow2 * r3pow3,
115 r2 * (3 * g3 + 6 * g4 * r3pow2 + g5 * r3pow4),
116 15 * g3 * r3 + 10 * g4 * r3pow3 + g5 * r3pow5};
117
118 auto D4 = SymTensor3d_4<T>{
119 3 * g2 + 6 * g3 * r1pow2 + g4 * r1pow4,
120 3 * g3 * r1 * r2 + g4 * r1pow3 * r2,
121 3 * g3 * r1 * r3 + g4 * r1pow3 * r3,
122 g2 + g3 * r1pow2 + g3 * r2pow2 + g4 * r1pow2 * r2pow2,
123 (g3 + g4 * r1pow2) * r2 * r3,
124 g2 + g3 * r1pow2 + g3 * r3pow2 + g4 * r1pow2 * r3pow2,
125 3 * g3 * r1 * r2 + g4 * r1 * r2pow3,
126 r1 * (g3 + g4 * r2pow2) * r3,
127 r1 * r2 * (g3 + g4 * r3pow2),
128 3 * g3 * r1 * r3 + g4 * r1 * r3pow3,
129 3 * g2 + 6 * g3 * r2pow2 + g4 * r2pow4,
130 3 * g3 * r2 * r3 + g4 * r2pow3 * r3,
131 g2 + g3 * r2pow2 + g3 * r3pow2 + g4 * r2pow2 * r3pow2,
132 3 * g3 * r2 * r3 + g4 * r2 * r3pow3,
133 3 * g2 + 6 * g3 * r3pow2 + g4 * r3pow4};
134
135 auto D3 = SymTensor3d_3<T>{
136 3 * g2 * r1 + g3 * r1pow3,
137 g2 * r2 + g3 * r1pow2 * r2,
138 g2 * r3 + g3 * r1pow2 * r3,
139 g2 * r1 + g3 * r1 * r2pow2,
140 g3 * r1 * r2 * r3,
141 g2 * r1 + g3 * r1 * r3pow2,
142 3 * g2 * r2 + g3 * r2pow3,
143 g2 * r3 + g3 * r2pow2 * r3,
144 g2 * r2 + g3 * r2 * r3pow2,
145 3 * g2 * r3 + g3 * r3pow3};
146
147 auto D2 = SymTensor3d_2<T>{
148 g1 + g2 * r1pow2,
149 g2 * r1 * r2,
150 g2 * r1 * r3,
151 g1 + g2 * r2pow2,
152 g2 * r2 * r3,
153 g1 + g2 * r3pow2};
154
155 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
156
157 auto D0 = g0;
158
159 return SymTensorCollection<T, 0, 5>{D0, D1, D2, D3, D4, D5};
160 }
161 };
162
163 template<class T>
164 class GreenFuncGravCartesian<T, 0, 4> {
165 public:
166 inline static shammath::SymTensorCollection<T, 0, 4> get_der_tensors(
167 const sycl::vec<T, 3> &r) {
168 using namespace shammath;
169 T r1 = r.x();
170 T r2 = r.y();
171 T r3 = r.z();
172
173 T r1pow2 = r.x() * r.x();
174 T r2pow2 = r.y() * r.y();
175 T r3pow2 = r.z() * r.z();
176
177 T r1pow3 = r.x() * r1pow2;
178 T r2pow3 = r.y() * r2pow2;
179 T r3pow3 = r.z() * r3pow2;
180
181 T r1pow4 = r.x() * r1pow3;
182 T r2pow4 = r.y() * r2pow3;
183 T r3pow4 = r.z() * r3pow3;
184
185 T rsq = r1pow2 + r2pow2 + r3pow2;
186
187 T rnorm = sycl::sqrt(rsq);
188
189 T rm2 = 1 / (rsq);
190
191 T g0 = 1 / rnorm;
192 T g1 = -1 * rm2 * g0;
193 T g2 = -3 * rm2 * g1;
194 T g3 = -5 * rm2 * g2;
195 T g4 = -7 * rm2 * g3;
196
197 auto D4 = SymTensor3d_4<T>{
198 3 * g2 + 6 * g3 * r1pow2 + g4 * r1pow4,
199 3 * g3 * r1 * r2 + g4 * r1pow3 * r2,
200 3 * g3 * r1 * r3 + g4 * r1pow3 * r3,
201 g2 + g3 * r1pow2 + g3 * r2pow2 + g4 * r1pow2 * r2pow2,
202 (g3 + g4 * r1pow2) * r2 * r3,
203 g2 + g3 * r1pow2 + g3 * r3pow2 + g4 * r1pow2 * r3pow2,
204 3 * g3 * r1 * r2 + g4 * r1 * r2pow3,
205 r1 * (g3 + g4 * r2pow2) * r3,
206 r1 * r2 * (g3 + g4 * r3pow2),
207 3 * g3 * r1 * r3 + g4 * r1 * r3pow3,
208 3 * g2 + 6 * g3 * r2pow2 + g4 * r2pow4,
209 3 * g3 * r2 * r3 + g4 * r2pow3 * r3,
210 g2 + g3 * r2pow2 + g3 * r3pow2 + g4 * r2pow2 * r3pow2,
211 3 * g3 * r2 * r3 + g4 * r2 * r3pow3,
212 3 * g2 + 6 * g3 * r3pow2 + g4 * r3pow4};
213
214 auto D3 = SymTensor3d_3<T>{
215 3 * g2 * r1 + g3 * r1pow3,
216 g2 * r2 + g3 * r1pow2 * r2,
217 g2 * r3 + g3 * r1pow2 * r3,
218 g2 * r1 + g3 * r1 * r2pow2,
219 g3 * r1 * r2 * r3,
220 g2 * r1 + g3 * r1 * r3pow2,
221 3 * g2 * r2 + g3 * r2pow3,
222 g2 * r3 + g3 * r2pow2 * r3,
223 g2 * r2 + g3 * r2 * r3pow2,
224 3 * g2 * r3 + g3 * r3pow3};
225
226 auto D2 = SymTensor3d_2<T>{
227 g1 + g2 * r1pow2,
228 g2 * r1 * r2,
229 g2 * r1 * r3,
230 g1 + g2 * r2pow2,
231 g2 * r2 * r3,
232 g1 + g2 * r3pow2};
233
234 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
235
236 auto D0 = g0;
237
238 return SymTensorCollection<T, 0, 4>{D0, D1, D2, D3, D4};
239 }
240 };
241
242 template<class T>
243 class GreenFuncGravCartesian<T, 0, 3> {
244 public:
245 inline static shammath::SymTensorCollection<T, 0, 3> get_der_tensors(
246 const sycl::vec<T, 3> &r) {
247 using namespace shammath;
248
249 T r1 = r.x();
250 T r2 = r.y();
251 T r3 = r.z();
252
253 T r1pow2 = r.x() * r.x();
254 T r2pow2 = r.y() * r.y();
255 T r3pow2 = r.z() * r.z();
256
257 T r1pow3 = r.x() * r1pow2;
258 T r2pow3 = r.y() * r2pow2;
259 T r3pow3 = r.z() * r3pow2;
260
261 T rsq = r1pow2 + r2pow2 + r3pow2;
262
263 T rnorm = sycl::sqrt(rsq);
264
265 T rm2 = 1 / (rsq);
266
267 T g0 = 1 / rnorm;
268 T g1 = -1 * rm2 * g0;
269 T g2 = -3 * rm2 * g1;
270 T g3 = -5 * rm2 * g2;
271
272 auto D3 = SymTensor3d_3<T>{
273 3 * g2 * r1 + g3 * r1pow3,
274 g2 * r2 + g3 * r1pow2 * r2,
275 g2 * r3 + g3 * r1pow2 * r3,
276 g2 * r1 + g3 * r1 * r2pow2,
277 g3 * r1 * r2 * r3,
278 g2 * r1 + g3 * r1 * r3pow2,
279 3 * g2 * r2 + g3 * r2pow3,
280 g2 * r3 + g3 * r2pow2 * r3,
281 g2 * r2 + g3 * r2 * r3pow2,
282 3 * g2 * r3 + g3 * r3pow3};
283
284 auto D2 = SymTensor3d_2<T>{
285 g1 + g2 * r1pow2,
286 g2 * r1 * r2,
287 g2 * r1 * r3,
288 g1 + g2 * r2pow2,
289 g2 * r2 * r3,
290 g1 + g2 * r3pow2};
291
292 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
293
294 auto D0 = g0;
295
296 return SymTensorCollection<T, 0, 3>{D0, D1, D2, D3};
297 }
298 };
299
300 template<class T>
301 class GreenFuncGravCartesian<T, 0, 2> {
302 public:
303 inline static shammath::SymTensorCollection<T, 0, 2> get_der_tensors(
304 const sycl::vec<T, 3> &r) {
305
306 using namespace shammath;
307
308 T r1 = r.x();
309 T r2 = r.y();
310 T r3 = r.z();
311
312 T r1pow2 = r.x() * r.x();
313 T r2pow2 = r.y() * r.y();
314 T r3pow2 = r.z() * r.z();
315
316 T r1pow3 = r.x() * r1pow2;
317 T r2pow3 = r.y() * r2pow2;
318 T r3pow3 = r.z() * r3pow2;
319
320 T rsq = r1pow2 + r2pow2 + r3pow2;
321
322 T rnorm = sycl::sqrt(rsq);
323
324 T rm2 = 1 / (rsq);
325
326 T g0 = 1 / rnorm;
327 T g1 = -1 * rm2 * g0;
328 T g2 = -3 * rm2 * g1;
329
330 auto D2 = SymTensor3d_2<T>{
331 g1 + g2 * r1pow2,
332 g2 * r1 * r2,
333 g2 * r1 * r3,
334 g1 + g2 * r2pow2,
335 g2 * r2 * r3,
336 g1 + g2 * r3pow2};
337
338 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
339
340 auto D0 = g0;
341
342 return SymTensorCollection<T, 0, 2>{D0, D1, D2};
343 }
344 };
345
346 template<class T>
347 class GreenFuncGravCartesian<T, 0, 1> {
348 public:
349 inline static shammath::SymTensorCollection<T, 0, 1> get_der_tensors(
350 const sycl::vec<T, 3> &r) {
351 using namespace shammath;
352 T r1 = r.x();
353 T r2 = r.y();
354 T r3 = r.z();
355
356 T r1pow2 = r.x() * r.x();
357 T r2pow2 = r.y() * r.y();
358 T r3pow2 = r.z() * r.z();
359
360 T r1pow3 = r.x() * r1pow2;
361 T r2pow3 = r.y() * r2pow2;
362 T r3pow3 = r.z() * r3pow2;
363
364 T rsq = r1pow2 + r2pow2 + r3pow2;
365
366 T rnorm = sycl::sqrt(rsq);
367
368 T rm2 = 1 / (rsq);
369
370 T g0 = 1 / rnorm;
371 T g1 = -1 * rm2 * g0;
372
373 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
374
375 auto D0 = g0;
376
377 return SymTensorCollection<T, 0, 1>{D0, D1};
378 }
379 };
380
381 template<class T>
382 class GreenFuncGravCartesian<T, 0, 0> {
383 public:
384 inline static shammath::SymTensorCollection<T, 0, 0> get_der_tensors(
385 const sycl::vec<T, 3> &r) {
386 using namespace shammath;
387 T r1 = r.x();
388 T r2 = r.y();
389 T r3 = r.z();
390
391 T r1pow2 = r.x() * r.x();
392 T r2pow2 = r.y() * r.y();
393 T r3pow2 = r.z() * r.z();
394
395 T r1pow3 = r.x() * r1pow2;
396 T r2pow3 = r.y() * r2pow2;
397 T r3pow3 = r.z() * r3pow2;
398
399 T rsq = r1pow2 + r2pow2 + r3pow2;
400
401 T rnorm = sycl::sqrt(rsq);
402
403 T rm2 = 1 / (rsq);
404
405 T g0 = 1 / rnorm;
406 T g1 = -1 * rm2 * g0;
407
408 auto D0 = g0;
409
410 return SymTensorCollection<T, 0, 0>{D0};
411 }
412 };
413
414 template<class T>
415 class GreenFuncGravCartesian<T, 1, 5> {
416 public:
417 inline static shammath::SymTensorCollection<T, 1, 5> get_der_tensors(
418 const sycl::vec<T, 3> &r) {
419 using namespace shammath;
420 T r1 = r.x();
421 T r2 = r.y();
422 T r3 = r.z();
423
424 T r1pow2 = r.x() * r.x();
425 T r2pow2 = r.y() * r.y();
426 T r3pow2 = r.z() * r.z();
427
428 T r1pow3 = r.x() * r1pow2;
429 T r2pow3 = r.y() * r2pow2;
430 T r3pow3 = r.z() * r3pow2;
431
432 T r1pow4 = r.x() * r1pow3;
433 T r2pow4 = r.y() * r2pow3;
434 T r3pow4 = r.z() * r3pow3;
435
436 T r1pow5 = r.x() * r1pow4;
437 T r2pow5 = r.y() * r2pow4;
438 T r3pow5 = r.z() * r3pow4;
439
440 T rsq = r1pow2 + r2pow2 + r3pow2;
441
442 T rnorm = sycl::sqrt(rsq);
443
444 T rm2 = 1 / (rsq);
445
446 T g0 = 1 / rnorm;
447 T g1 = -1 * rm2 * g0;
448 T g2 = -3 * rm2 * g1;
449 T g3 = -5 * rm2 * g2;
450 T g4 = -7 * rm2 * g3;
451 T g5 = -9 * rm2 * g4;
452
453 auto D5 = SymTensor3d_5<T>{
454 15 * g3 * r1 + 10 * g4 * r1pow3 + g5 * r1pow5,
455 (3 * g3 + 6 * g4 * r1pow2 + g5 * r1pow4) * r2,
456 (3 * g3 + 6 * g4 * r1pow2 + g5 * r1pow4) * r3,
457 3 * g3 * r1 + g4 * r1pow3 + 3 * g4 * r1 * r2pow2 + g5 * r1pow3 * r2pow2,
458 (3 * g4 * r1 + g5 * r1pow3) * r2 * r3,
459 3 * g3 * r1 + g4 * r1pow3 + 3 * g4 * r1 * r3pow2 + g5 * r1pow3 * r3pow2,
460 3 * g3 * r2 + 3 * g4 * r1pow2 * r2 + g4 * r2pow3 + g5 * r1pow2 * r2pow3,
461 (g3 + g4 * r1pow2 + g4 * r2pow2 + g5 * r1pow2 * r2pow2) * r3,
462 r2 * (g3 + g4 * r1pow2 + g4 * r3pow2 + g5 * r1pow2 * r3pow2),
463 3 * g3 * r3 + 3 * g4 * r1pow2 * r3 + g4 * r3pow3 + g5 * r1pow2 * r3pow3,
464 r1 * (3 * g3 + 6 * g4 * r2pow2 + g5 * r2pow4),
465 r1 * (3 * g4 * r2 + g5 * r2pow3) * r3,
466 r1 * (g3 + g4 * r2pow2 + g4 * r3pow2 + g5 * r2pow2 * r3pow2),
467 r1 * r2 * (3 * g4 * r3 + g5 * r3pow3),
468 r1 * (3 * g3 + 6 * g4 * r3pow2 + g5 * r3pow4),
469 15 * g3 * r2 + 10 * g4 * r2pow3 + g5 * r2pow5,
470 (3 * g3 + 6 * g4 * r2pow2 + g5 * r2pow4) * r3,
471 3 * g3 * r2 + g4 * r2pow3 + 3 * g4 * r2 * r3pow2 + g5 * r2pow3 * r3pow2,
472 3 * g3 * r3 + 3 * g4 * r2pow2 * r3 + g4 * r3pow3 + g5 * r2pow2 * r3pow3,
473 r2 * (3 * g3 + 6 * g4 * r3pow2 + g5 * r3pow4),
474 15 * g3 * r3 + 10 * g4 * r3pow3 + g5 * r3pow5};
475
476 auto D4 = SymTensor3d_4<T>{
477 3 * g2 + 6 * g3 * r1pow2 + g4 * r1pow4,
478 3 * g3 * r1 * r2 + g4 * r1pow3 * r2,
479 3 * g3 * r1 * r3 + g4 * r1pow3 * r3,
480 g2 + g3 * r1pow2 + g3 * r2pow2 + g4 * r1pow2 * r2pow2,
481 (g3 + g4 * r1pow2) * r2 * r3,
482 g2 + g3 * r1pow2 + g3 * r3pow2 + g4 * r1pow2 * r3pow2,
483 3 * g3 * r1 * r2 + g4 * r1 * r2pow3,
484 r1 * (g3 + g4 * r2pow2) * r3,
485 r1 * r2 * (g3 + g4 * r3pow2),
486 3 * g3 * r1 * r3 + g4 * r1 * r3pow3,
487 3 * g2 + 6 * g3 * r2pow2 + g4 * r2pow4,
488 3 * g3 * r2 * r3 + g4 * r2pow3 * r3,
489 g2 + g3 * r2pow2 + g3 * r3pow2 + g4 * r2pow2 * r3pow2,
490 3 * g3 * r2 * r3 + g4 * r2 * r3pow3,
491 3 * g2 + 6 * g3 * r3pow2 + g4 * r3pow4};
492
493 auto D3 = SymTensor3d_3<T>{
494 3 * g2 * r1 + g3 * r1pow3,
495 g2 * r2 + g3 * r1pow2 * r2,
496 g2 * r3 + g3 * r1pow2 * r3,
497 g2 * r1 + g3 * r1 * r2pow2,
498 g3 * r1 * r2 * r3,
499 g2 * r1 + g3 * r1 * r3pow2,
500 3 * g2 * r2 + g3 * r2pow3,
501 g2 * r3 + g3 * r2pow2 * r3,
502 g2 * r2 + g3 * r2 * r3pow2,
503 3 * g2 * r3 + g3 * r3pow3};
504
505 auto D2 = SymTensor3d_2<T>{
506 g1 + g2 * r1pow2,
507 g2 * r1 * r2,
508 g2 * r1 * r3,
509 g1 + g2 * r2pow2,
510 g2 * r2 * r3,
511 g1 + g2 * r3pow2};
512
513 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
514
515 return SymTensorCollection<T, 1, 5>{D1, D2, D3, D4, D5};
516 }
517 };
518
519 template<class T>
520 class GreenFuncGravCartesian<T, 1, 4> {
521 public:
522 inline static shammath::SymTensorCollection<T, 1, 4> get_der_tensors(
523 const sycl::vec<T, 3> &r) {
524
525 using namespace shammath;
526
527 T r1 = r.x();
528 T r2 = r.y();
529 T r3 = r.z();
530
531 T r1pow2 = r.x() * r.x();
532 T r2pow2 = r.y() * r.y();
533 T r3pow2 = r.z() * r.z();
534
535 T r1pow3 = r.x() * r1pow2;
536 T r2pow3 = r.y() * r2pow2;
537 T r3pow3 = r.z() * r3pow2;
538
539 T r1pow4 = r.x() * r1pow3;
540 T r2pow4 = r.y() * r2pow3;
541 T r3pow4 = r.z() * r3pow3;
542
543 T rsq = r1pow2 + r2pow2 + r3pow2;
544
545 T rnorm = sycl::sqrt(rsq);
546
547 T rm2 = 1 / (rsq);
548
549 T g0 = 1 / rnorm;
550 T g1 = -1 * rm2 * g0;
551 T g2 = -3 * rm2 * g1;
552 T g3 = -5 * rm2 * g2;
553 T g4 = -7 * rm2 * g3;
554
555 auto D4 = SymTensor3d_4<T>{
556 3 * g2 + 6 * g3 * r1pow2 + g4 * r1pow4,
557 3 * g3 * r1 * r2 + g4 * r1pow3 * r2,
558 3 * g3 * r1 * r3 + g4 * r1pow3 * r3,
559 g2 + g3 * r1pow2 + g3 * r2pow2 + g4 * r1pow2 * r2pow2,
560 (g3 + g4 * r1pow2) * r2 * r3,
561 g2 + g3 * r1pow2 + g3 * r3pow2 + g4 * r1pow2 * r3pow2,
562 3 * g3 * r1 * r2 + g4 * r1 * r2pow3,
563 r1 * (g3 + g4 * r2pow2) * r3,
564 r1 * r2 * (g3 + g4 * r3pow2),
565 3 * g3 * r1 * r3 + g4 * r1 * r3pow3,
566 3 * g2 + 6 * g3 * r2pow2 + g4 * r2pow4,
567 3 * g3 * r2 * r3 + g4 * r2pow3 * r3,
568 g2 + g3 * r2pow2 + g3 * r3pow2 + g4 * r2pow2 * r3pow2,
569 3 * g3 * r2 * r3 + g4 * r2 * r3pow3,
570 3 * g2 + 6 * g3 * r3pow2 + g4 * r3pow4};
571
572 auto D3 = SymTensor3d_3<T>{
573 3 * g2 * r1 + g3 * r1pow3,
574 g2 * r2 + g3 * r1pow2 * r2,
575 g2 * r3 + g3 * r1pow2 * r3,
576 g2 * r1 + g3 * r1 * r2pow2,
577 g3 * r1 * r2 * r3,
578 g2 * r1 + g3 * r1 * r3pow2,
579 3 * g2 * r2 + g3 * r2pow3,
580 g2 * r3 + g3 * r2pow2 * r3,
581 g2 * r2 + g3 * r2 * r3pow2,
582 3 * g2 * r3 + g3 * r3pow3};
583
584 auto D2 = SymTensor3d_2<T>{
585 g1 + g2 * r1pow2,
586 g2 * r1 * r2,
587 g2 * r1 * r3,
588 g1 + g2 * r2pow2,
589 g2 * r2 * r3,
590 g1 + g2 * r3pow2};
591
592 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
593
594 return SymTensorCollection<T, 1, 4>{D1, D2, D3, D4};
595 }
596 };
597
598 template<class T>
599 class GreenFuncGravCartesian<T, 1, 3> {
600 public:
601 inline static shammath::SymTensorCollection<T, 1, 3> get_der_tensors(
602 const sycl::vec<T, 3> &r) {
603
604 using namespace shammath;
605
606 T r1 = r.x();
607 T r2 = r.y();
608 T r3 = r.z();
609
610 T r1pow2 = r.x() * r.x();
611 T r2pow2 = r.y() * r.y();
612 T r3pow2 = r.z() * r.z();
613
614 T r1pow3 = r.x() * r1pow2;
615 T r2pow3 = r.y() * r2pow2;
616 T r3pow3 = r.z() * r3pow2;
617
618 T rsq = r1pow2 + r2pow2 + r3pow2;
619
620 T rnorm = sycl::sqrt(rsq);
621
622 T rm2 = 1 / (rsq);
623
624 T g0 = 1 / rnorm;
625 T g1 = -1 * rm2 * g0;
626 T g2 = -3 * rm2 * g1;
627 T g3 = -5 * rm2 * g2;
628
629 auto D3 = SymTensor3d_3<T>{
630 3 * g2 * r1 + g3 * r1pow3,
631 g2 * r2 + g3 * r1pow2 * r2,
632 g2 * r3 + g3 * r1pow2 * r3,
633 g2 * r1 + g3 * r1 * r2pow2,
634 g3 * r1 * r2 * r3,
635 g2 * r1 + g3 * r1 * r3pow2,
636 3 * g2 * r2 + g3 * r2pow3,
637 g2 * r3 + g3 * r2pow2 * r3,
638 g2 * r2 + g3 * r2 * r3pow2,
639 3 * g2 * r3 + g3 * r3pow3};
640
641 auto D2 = SymTensor3d_2<T>{
642 g1 + g2 * r1pow2,
643 g2 * r1 * r2,
644 g2 * r1 * r3,
645 g1 + g2 * r2pow2,
646 g2 * r2 * r3,
647 g1 + g2 * r3pow2};
648
649 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
650
651 return SymTensorCollection<T, 1, 3>{D1, D2, D3};
652 }
653 };
654
655 template<class T>
656 class GreenFuncGravCartesian<T, 1, 2> {
657 public:
658 inline static shammath::SymTensorCollection<T, 1, 2> get_der_tensors(
659 const sycl::vec<T, 3> &r) {
660
661 using namespace shammath;
662
663 T r1 = r.x();
664 T r2 = r.y();
665 T r3 = r.z();
666
667 T r1pow2 = r.x() * r.x();
668 T r2pow2 = r.y() * r.y();
669 T r3pow2 = r.z() * r.z();
670
671 T r1pow3 = r.x() * r1pow2;
672 T r2pow3 = r.y() * r2pow2;
673 T r3pow3 = r.z() * r3pow2;
674
675 T rsq = r1pow2 + r2pow2 + r3pow2;
676
677 T rnorm = sycl::sqrt(rsq);
678
679 T rm2 = 1 / (rsq);
680
681 T g0 = 1 / rnorm;
682 T g1 = -1 * rm2 * g0;
683 T g2 = -3 * rm2 * g1;
684
685 auto D2 = SymTensor3d_2<T>{
686 g1 + g2 * r1pow2,
687 g2 * r1 * r2,
688 g2 * r1 * r3,
689 g1 + g2 * r2pow2,
690 g2 * r2 * r3,
691 g1 + g2 * r3pow2};
692
693 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
694
695 return SymTensorCollection<T, 1, 2>{D1, D2};
696 }
697 };
698
699 template<class T>
700 class GreenFuncGravCartesian<T, 1, 1> {
701 public:
702 inline static shammath::SymTensorCollection<T, 1, 1> get_der_tensors(
703 const sycl::vec<T, 3> &r) {
704 using namespace shammath;
705 T r1 = r.x();
706 T r2 = r.y();
707 T r3 = r.z();
708
709 T r1pow2 = r.x() * r.x();
710 T r2pow2 = r.y() * r.y();
711 T r3pow2 = r.z() * r.z();
712
713 T r1pow3 = r.x() * r1pow2;
714 T r2pow3 = r.y() * r2pow2;
715 T r3pow3 = r.z() * r3pow2;
716
717 T rsq = r1pow2 + r2pow2 + r3pow2;
718
719 T rnorm = sycl::sqrt(rsq);
720
721 T rm2 = 1 / (rsq);
722
723 T g0 = 1 / rnorm;
724 T g1 = -1 * rm2 * g0;
725
726 auto D1 = SymTensor3d_1<T>{g1 * r1, g1 * r2, g1 * r3};
727
728 return SymTensorCollection<T, 1, 1>{D1};
729 }
730 };
731
732#endif
733} // namespace shamphys
aliases_int.hpp
shamphys::GreenFuncGravCartesian
Utility to get the derivatives of the Green function for gravity in Cartesian coordinates.
Definition GreenFuncGravCartesian.hpp:37
shammath
namespace for math utility
Definition AABB.hpp:26
shammath::SymTensor3d_1
Definition symtensors.hpp:23
shammath::SymTensor3d_2
Definition symtensors.hpp:91
shammath::SymTensor3d_3
Definition symtensors.hpp:192
shammath::SymTensor3d_4
Definition symtensors.hpp:354
shammath::SymTensor3d_5
Definition symtensors.hpp:589
symtensor_collections.hpp
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8