Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
grav_moment_offset.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
17
18#include "shambase/aliases_int.hpp"
19#include "shambase/type_traits.hpp"
20#include "shammath/symtensor_collections.hpp"
21#include "shamphys/fmm/grav_moments.hpp"
22
23namespace shamphys {
24
25 // Offset the gravitational moment derivatives
26 template<class T, u32 low_order, u32 high_order>
27 inline shammath::SymTensorCollection<T, low_order, high_order> offset_dM_mat_delta(
28 const shammath::SymTensorCollection<T, low_order, high_order> &dM,
29 const sycl::vec<T, 3> &offset) {
30
31 using namespace shammath;
32
33 static constexpr T inv_factorial_0 = 1. / 1;
34 static constexpr T inv_factorial_1 = 1. / 1;
35 static constexpr T inv_factorial_2 = 1. / 2;
36 static constexpr T inv_factorial_3 = 1. / 6;
37 static constexpr T inv_factorial_4 = 1. / 24;
38 static constexpr T inv_factorial_5 = 1. / 120;
39
40 if constexpr (low_order == 0 && high_order == 5) {
41
42 SymTensorCollection<T, 0, 5> h = SymTensorCollection<T, 0, 5>::from_vec(offset);
43
44 auto &h_0 = h.t0;
45 auto &h_1 = h.t1;
46 auto &h_2 = h.t2;
47 auto &h_3 = h.t3;
48 auto &h_4 = h.t4;
49 auto &h_5 = h.t5;
50
51 auto &dM_0 = dM.t0;
52 auto &dM_1 = dM.t1;
53 auto &dM_2 = dM.t2;
54 auto &dM_3 = dM.t3;
55 auto &dM_4 = dM.t4;
56 auto &dM_5 = dM.t5;
57
58 shammath::SymTensorCollection<T, 0, 5> dM_ret;
59
60 auto &dM_ret_0 = dM_ret.t0;
61 auto &dM_ret_1 = dM_ret.t1;
62 auto &dM_ret_2 = dM_ret.t2;
63 auto &dM_ret_3 = dM_ret.t3;
64 auto &dM_ret_4 = dM_ret.t4;
65 auto &dM_ret_5 = dM_ret.t5;
66
67 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
68
69 dM_ret_0 = inv_factorial_0 * h_0 * dM_0 + inv_factorial_1 * h_1 * dM_1
70 + inv_factorial_2 * h_2 * dM_2 + inv_factorial_3 * h_3 * dM_3
71 + inv_factorial_4 * h_4 * dM_4 + inv_factorial_5 * h_5 * dM_5;
72
73 dM_ret_1 = inv_factorial_0 * h_0 * dM_1 + inv_factorial_1 * h_1 * dM_2
74 + inv_factorial_2 * h_2 * dM_3 + inv_factorial_3 * h_3 * dM_4
75 + inv_factorial_4 * h_4 * dM_5;
76
77 dM_ret_2 = inv_factorial_0 * h_0 * dM_2 + inv_factorial_1 * h_1 * dM_3
78 + inv_factorial_2 * h_2 * dM_4 + inv_factorial_3 * h_3 * dM_5;
79
80 dM_ret_3 = inv_factorial_0 * h_0 * dM_3 + inv_factorial_1 * h_1 * dM_4
81 + inv_factorial_2 * h_2 * dM_5;
82
83 dM_ret_4 = inv_factorial_0 * h_0 * dM_4 + inv_factorial_1 * h_1 * dM_5;
84 dM_ret_5 = inv_factorial_0 * h_0 * dM_5;
85
86 return dM_ret;
87 } else if constexpr (low_order == 1 && high_order == 5) {
88
89 SymTensorCollection<T, 0, 5> h = SymTensorCollection<T, 0, 5>::from_vec(offset);
90
91 auto &h_0 = h.t0;
92 auto &h_1 = h.t1;
93 auto &h_2 = h.t2;
94 auto &h_3 = h.t3;
95 auto &h_4 = h.t4;
96 auto &h_5 = h.t5;
97
98 auto &dM_1 = dM.t1;
99 auto &dM_2 = dM.t2;
100 auto &dM_3 = dM.t3;
101 auto &dM_4 = dM.t4;
102 auto &dM_5 = dM.t5;
103
104 shammath::SymTensorCollection<T, 1, 5> dM_ret;
105
106 auto &dM_ret_1 = dM_ret.t1;
107 auto &dM_ret_2 = dM_ret.t2;
108 auto &dM_ret_3 = dM_ret.t3;
109 auto &dM_ret_4 = dM_ret.t4;
110 auto &dM_ret_5 = dM_ret.t5;
111
112 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
113
114 dM_ret_1 = inv_factorial_0 * h_0 * dM_1 + inv_factorial_1 * h_1 * dM_2
115 + inv_factorial_2 * h_2 * dM_3 + inv_factorial_3 * h_3 * dM_4
116 + inv_factorial_4 * h_4 * dM_5;
117
118 dM_ret_2 = inv_factorial_0 * h_0 * dM_2 + inv_factorial_1 * h_1 * dM_3
119 + inv_factorial_2 * h_2 * dM_4 + inv_factorial_3 * h_3 * dM_5;
120
121 dM_ret_3 = inv_factorial_0 * h_0 * dM_3 + inv_factorial_1 * h_1 * dM_4
122 + inv_factorial_2 * h_2 * dM_5;
123
124 dM_ret_4 = inv_factorial_0 * h_0 * dM_4 + inv_factorial_1 * h_1 * dM_5;
125 dM_ret_5 = inv_factorial_0 * h_0 * dM_5;
126
127 return dM_ret;
128 } else if constexpr (low_order == 1 && high_order == 4) {
129
130 SymTensorCollection<T, 0, 4> h = SymTensorCollection<T, 0, 4>::from_vec(offset);
131
132 auto &h_0 = h.t0;
133 auto &h_1 = h.t1;
134 auto &h_2 = h.t2;
135 auto &h_3 = h.t3;
136 auto &h_4 = h.t4;
137
138 auto &dM_1 = dM.t1;
139 auto &dM_2 = dM.t2;
140 auto &dM_3 = dM.t3;
141 auto &dM_4 = dM.t4;
142
143 shammath::SymTensorCollection<T, 1, 4> dM_ret;
144
145 auto &dM_ret_1 = dM_ret.t1;
146 auto &dM_ret_2 = dM_ret.t2;
147 auto &dM_ret_3 = dM_ret.t3;
148 auto &dM_ret_4 = dM_ret.t4;
149
150 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
151
152 dM_ret_1 = inv_factorial_0 * h_0 * dM_1 + inv_factorial_1 * h_1 * dM_2
153 + inv_factorial_2 * h_2 * dM_3 + inv_factorial_3 * h_3 * dM_4;
154
155 dM_ret_2 = inv_factorial_0 * h_0 * dM_2 + inv_factorial_1 * h_1 * dM_3
156 + inv_factorial_2 * h_2 * dM_4;
157
158 dM_ret_3 = inv_factorial_0 * h_0 * dM_3 + inv_factorial_1 * h_1 * dM_4;
159
160 dM_ret_4 = inv_factorial_0 * h_0 * dM_4;
161
162 return dM_ret;
163 } else if constexpr (low_order == 1 && high_order == 3) {
164
165 SymTensorCollection<T, 0, 3> h = SymTensorCollection<T, 0, 3>::from_vec(offset);
166
167 auto &h_0 = h.t0;
168 auto &h_1 = h.t1;
169 auto &h_2 = h.t2;
170 auto &h_3 = h.t3;
171
172 auto &dM_1 = dM.t1;
173 auto &dM_2 = dM.t2;
174 auto &dM_3 = dM.t3;
175
176 shammath::SymTensorCollection<T, 1, 3> dM_ret;
177
178 auto &dM_ret_1 = dM_ret.t1;
179 auto &dM_ret_2 = dM_ret.t2;
180 auto &dM_ret_3 = dM_ret.t3;
181
182 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
183
184 dM_ret_1 = inv_factorial_0 * h_0 * dM_1 + inv_factorial_1 * h_1 * dM_2
185 + inv_factorial_2 * h_2 * dM_3;
186
187 dM_ret_2 = inv_factorial_0 * h_0 * dM_2 + inv_factorial_1 * h_1 * dM_3;
188
189 dM_ret_3 = inv_factorial_0 * h_0 * dM_3;
190
191 return dM_ret;
192 } else if constexpr (low_order == 1 && high_order == 2) {
193
194 SymTensorCollection<T, 0, 1> h = SymTensorCollection<T, 0, 1>::from_vec(offset);
195
196 auto &h_0 = h.t0;
197 auto &h_1 = h.t1;
198
199 auto &dM_1 = dM.t1;
200 auto &dM_2 = dM.t2;
201
202 shammath::SymTensorCollection<T, 1, 2> dM_ret;
203
204 auto &dM_ret_1 = dM_ret.t1;
205 auto &dM_ret_2 = dM_ret.t2;
206
207 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
208
209 dM_ret_1 = inv_factorial_0 * h_0 * dM_1 + inv_factorial_1 * h_1 * dM_2;
210
211 dM_ret_2 = inv_factorial_0 * h_0 * dM_2;
212
213 return dM_ret;
214 } else if constexpr (low_order == 1 && high_order == 1) {
215 T h_0 = 1;
216
217 auto &dM_1 = dM.t1;
218
219 shammath::SymTensorCollection<T, 1, 1> dM_ret;
220
221 auto &dM_ret_1 = dM_ret.t1;
222
223 // dM_k = sum_{l=k}^p \frac{1}{(l-k)!} h^(l-k).dM_l
224
225 dM_ret_1 = inv_factorial_0 * h_0 * dM_1;
226
227 return dM_ret;
228 } else {
229 static_assert(shambase::always_false_v<T>, "This combination of orders is not valid");
230 }
231 }
232
233 // Offset the gravitational moment derivatives (with two vecs)
234 template<class T, u32 low_order, u32 high_order>
235 inline shammath::SymTensorCollection<T, low_order, high_order> offset_dM_mat(
236 const shammath::SymTensorCollection<T, low_order, high_order> &dM,
237 const sycl::vec<T, 3> &from,
238 const sycl::vec<T, 3> &to) {
239 return offset_dM_mat_delta(dM, to - from);
240 }
241
242} // namespace shamphys
aliases_int.hpp
grav_moments.hpp
shambase::always_false_v
constexpr bool always_false_v
Helper variable template that is always false. Especially useful to perform static asserts based on t...
Definition type_traits.hpp:234
type_traits.hpp
Traits for C++ types.
symtensor_collections.hpp
Generated on for Shamrock by doxygen 1.15.0