17 #ifndef SHAPEFACTORS_H_
18 #define SHAPEFACTORS_H_
27 template <
int depos_order>
44 using namespace amrex::literals;
46 const auto j =
static_cast<int>(amrex::Math::floor(xmid+0.5_rt));
60 using namespace amrex::literals;
62 const auto j =
static_cast<int>(amrex::Math::floor(xmid));
63 const amrex::Real xint = xmid-
j;
64 sx[0] = 1.0_rt - xint;
78 using namespace amrex::literals;
80 const auto j =
static_cast<int>(amrex::Math::floor(xmid+0.5_rt));
81 const amrex::Real xint = xmid-
j;
82 sx[0] = 0.5_rt*(0.5_rt-xint)*(0.5_rt-xint);
83 sx[1] = 0.75_rt-xint*xint;
84 sx[2] = 0.5_rt*(0.5_rt+xint)*(0.5_rt+xint);
98 using namespace amrex::literals;
100 const auto j =
static_cast<int>(amrex::Math::floor(xmid));
101 const amrex::Real xint = xmid-
j;
102 sx[0] = 1.0_rt/6.0_rt*(1.0_rt-xint)*(1.0_rt-xint)*(1.0_rt-xint);
103 sx[1] = 2.0_rt/3.0_rt-xint*xint*(1.0_rt-xint/2.0_rt);
104 sx[2] = 2.0_rt/3.0_rt-(1.0_rt-xint)*(1.0_rt-xint)*(1.0_rt-0.5_rt*(1.0_rt-xint));
105 sx[3] = 1.0_rt/6.0_rt*xint*xint*xint;
125 using namespace amrex::literals;
126 if constexpr (branchless) {
127 if constexpr (depos_order==0) {
128 return {1._rt,
static_cast<int>(amrex::Math::floor(xmid+0.5_rt))};
129 }
else if constexpr (depos_order==1) {
130 const amrex::Real xfloor = amrex::Math::floor(xmid);
131 const amrex::Real xint = xmid-xfloor;
132 const amrex::Real ic =
static_cast<amrex::Real
>(ix);
133 const amrex::Real shape = (1._rt-ic) + (2._rt*ic-1._rt) * xint;
134 return {shape,
static_cast<int>(xfloor)+ix};
135 }
else if constexpr (depos_order==2) {
136 const amrex::Real xfloor = amrex::Math::floor(xmid+0.5_rt);
137 const amrex::Real xint = xmid-xfloor;
138 const int ic = ix - 1;
139 const amrex::Real icr =
static_cast<amrex::Real
>(ic);
140 const amrex::Real shape = (-0.625_rt*icr*icr + 0.75_rt)
141 + (0.5_rt*icr + 1.5_rt*icr*icr*xint - xint) * xint;
142 return {shape,
static_cast<int>(xfloor)+ic};
143 }
else if constexpr (depos_order==3) {
144 const amrex::Real xfloor = amrex::Math::floor(xmid);
145 const amrex::Real xint = xmid-xfloor;
146 const int ic = ix + 1;
147 const amrex::Real icr =
static_cast<amrex::Real
>(ic&1);
148 const amrex::Real xint_s = icr + (1._rt - 2._rt*icr) * xint;
149 const amrex::Real cond =
static_cast<amrex::Real
>(ic&2);
150 const amrex::Real shape = (1.0_rt/6.0_rt)*(cond*2._rt
151 + xint_s*xint_s*(cond*(xint_s-3._rt) + xint_s));
152 return {shape,
static_cast<int>(xfloor)+ix-1};
155 if constexpr (depos_order==0) {
156 return {1._rt,
static_cast<int>(amrex::Math::floor(xmid+0.5_rt))};
157 }
else if constexpr (depos_order==1) {
158 const amrex::Real xfloor = amrex::Math::floor(xmid);
159 const amrex::Real xint = xmid-xfloor;
161 return {1.0_rt - xint,
static_cast<int>(xfloor)};
163 return {xint,
static_cast<int>(xfloor)+1};
165 }
else if constexpr (depos_order==2) {
166 const amrex::Real xfloor = amrex::Math::floor(xmid+0.5_rt);
167 const amrex::Real xint = xmid-xfloor;
169 return {0.5_rt*(0.5_rt-xint)*(0.5_rt-xint),
static_cast<int>(xfloor)-1};
171 return {0.75_rt-xint*xint,
static_cast<int>(xfloor)};
173 return {0.5_rt*(0.5_rt+xint)*(0.5_rt+xint),
static_cast<int>(xfloor)+1};
175 }
else if constexpr (depos_order==3) {
176 const amrex::Real xfloor = amrex::Math::floor(xmid);
177 const amrex::Real xint = xmid-xfloor;
179 return {1.0_rt/6.0_rt*(1.0_rt-xint)*(1.0_rt-xint)*(1.0_rt-xint),
180 static_cast<int>(xfloor)-1};
182 return {2.0_rt/3.0_rt-xint*xint*(1.0_rt-xint/2.0_rt),
183 static_cast<int>(xfloor)};
185 return {2.0_rt/3.0_rt-(1.0_rt-xint)*(1.0_rt-xint)*(1.0_rt-0.5_rt*(1.0_rt-xint)),
186 static_cast<int>(xfloor)+1};
188 return {1.0_rt/6.0_rt*xint*xint*xint,
189 static_cast<int>(xfloor)+2};
193 static_assert(0 <= depos_order && depos_order <= 3);
214 using namespace amrex::literals;
215 if constexpr (derivative_type==0) {
216 if constexpr (depos_order==0) {
218 const amrex::Real xfloor = amrex::Math::floor(xmid);
219 const amrex::Real s_x = 1;
220 const amrex::Real sdx = 0;
221 return {s_x, -sdx,
static_cast<int>(xfloor)};
222 }
else if constexpr (depos_order==1) {
223 const amrex::Real xfloor = amrex::Math::floor(xmid);
224 const amrex::Real xint = xmid-xfloor;
226 const amrex::Real s_x = 1 - xint;
227 const amrex::Real sdx = -1;
228 return {s_x, -sdx,
static_cast<int>(xfloor)+ix};
230 const amrex::Real s_x = xint;
231 const amrex::Real sdx = 1;
232 return {s_x, -sdx,
static_cast<int>(xfloor)+ix};
234 }
else if constexpr (depos_order==2) {
236 const amrex::Real xfloor = amrex::Math::floor(xmid);
237 const amrex::Real xint = xmid-xfloor;
238 const amrex::Real xint_2 = xint*xint;
240 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_2 - 1.0_rt*xint + 0.5_rt;
241 const amrex::Real sdx = xint - 1.0_rt;
242 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
244 const amrex::Real s_x = -xint_2 + 1.0_rt*xint + 0.5_rt;
245 const amrex::Real sdx = 1.0_rt - 2*xint;
246 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
248 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_2;
249 const amrex::Real sdx = xint;
250 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
252 }
else if constexpr (depos_order==3) {
253 const amrex::Real xfloor = amrex::Math::floor(xmid);
254 const amrex::Real xint = xmid-xfloor;
255 const amrex::Real xint_2 = xint*xint;
256 const amrex::Real xint_3 = xint_2*xint;
258 const amrex::Real s_x = -1.0_rt/6.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint + 1.0_rt/6.0_rt;
259 const amrex::Real sdx = -1.0_rt/2.0_rt*xint_2 + xint - 1.0_rt/2.0_rt;
260 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
262 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_3 - xint_2 + 2.0_rt/3.0_rt;
263 const amrex::Real sdx = (3.0_rt/2.0_rt)*xint_2 - 2*xint;
264 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
266 const amrex::Real s_x = -1.0_rt/2.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 + (1.0_rt/2.0_rt)*xint + 1.0_rt/6.0_rt;
267 const amrex::Real sdx = -3.0_rt/2.0_rt*xint_2 + xint + 1.0_rt/2.0_rt;
268 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
270 const amrex::Real s_x = (1.0_rt/6.0_rt)*xint_3;
271 const amrex::Real sdx = (1.0_rt/2.0_rt)*xint_2;
272 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
275 }
else if constexpr (derivative_type==1) {
276 if constexpr (depos_order==0) {
277 const amrex::Real xfloor = amrex::Math::floor(xmid);
278 const amrex::Real xint = xmid-xfloor;
280 const amrex::Real s_x = (xint < 0.5_rt) ? 1 : 0;
281 const amrex::Real sdx = -1;
282 return {s_x, -sdx,
static_cast<int>(xfloor)+ix};
284 const amrex::Real s_x = (xint < 0.5_rt) ? 0 : 1;
285 const amrex::Real sdx = 1;
286 return {s_x, -sdx,
static_cast<int>(xfloor)+ix};
288 }
else if constexpr (depos_order==1) {
290 const amrex::Real xfloor = amrex::Math::floor(xmid);
291 const amrex::Real xint = xmid-xfloor;
293 const amrex::Real s_x = (xint < 0.5_rt) ? 1.0_rt/2.0_rt - xint : 0;
294 const amrex::Real sdx = xint - 1;
295 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
297 const amrex::Real s_x = (xint < 0.5_rt) ? xint + 1.0_rt/2.0_rt : 3.0_rt/2.0_rt - xint;
298 const amrex::Real sdx = 1 - 2*xint;
299 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
301 const amrex::Real s_x = (xint < 0.5_rt) ? 0 : xint - 1.0_rt/2.0_rt;
302 const amrex::Real sdx = xint;
303 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
305 }
else if constexpr (depos_order==2) {
306 const amrex::Real xfloor = amrex::Math::floor(xmid);
307 const amrex::Real xint = xmid-xfloor;
308 const amrex::Real xint_2 = xint*xint;
310 const amrex::Real s_x = (xint < 0.5_rt) ?
311 (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint + 1.0_rt/8.0_rt : 0;
312 const amrex::Real sdx = -1.0_rt/2.0_rt*xint_2 + xint - 1.0_rt/2.0_rt;
313 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
315 const amrex::Real s_x = (xint < 0.5_rt) ?
316 3.0_rt/4.0_rt - xint_2 : (1.0_rt/2.0_rt)*xint_2 - 3.0_rt/2.0_rt*xint + 9.0_rt/8.0_rt;
317 const amrex::Real sdx = (3.0_rt/2.0_rt)*xint_2 - 2*xint;
318 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
320 const amrex::Real s_x = (xint < 0.5_rt) ?
321 (1.0_rt/2.0_rt)*xint_2 + (1.0_rt/2.0_rt)*xint + 1.0_rt/8.0_rt : -xint_2 + 2*xint - 1.0_rt/4.0_rt;
322 const amrex::Real sdx = -3.0_rt/2.0_rt*xint_2 + xint + 1.0_rt/2.0_rt;
323 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
325 const amrex::Real s_x = (xint < 0.5_rt) ?
326 0 : (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint + 1.0_rt/8.0_rt;
327 const amrex::Real sdx = (1.0_rt/2.0_rt)*xint_2;
328 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
330 }
else if constexpr (depos_order==3) {
332 const amrex::Real xfloor = amrex::Math::floor(xmid);
333 const amrex::Real xint = xmid-xfloor;
334 const amrex::Real xint_2 = xint*xint;
335 const amrex::Real xint_3 = xint_2*xint;
337 const amrex::Real s_x = (xint < 0.5_rt) ?
338 -1.0_rt/6.0_rt*xint_3 + (1.0_rt/4.0_rt)*xint_2 - 1.0_rt/8.0_rt*xint + 1.0_rt/48.0_rt :
340 const amrex::Real sdx = (1.0_rt/6.0_rt)*xint_3 - 1.0_rt/2.0_rt*xint_2 + (1.0_rt/2.0_rt)*xint - 1.0_rt/6.0_rt;
341 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
343 const amrex::Real s_x = (xint < 0.5_rt) ?
344 (1.0_rt/2.0_rt)*xint_3 - 1.0_rt/4.0_rt*xint_2 - 5.0_rt/8.0_rt*xint + 23.0_rt/48.0_rt :
345 -1.0_rt/6.0_rt*xint_3 + (3.0_rt/4.0_rt)*xint_2 - 9.0_rt/8.0_rt*xint + 9.0_rt/16.0_rt;
346 const amrex::Real sdx = -2.0_rt/3.0_rt*xint_3 + (3.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint - 1.0_rt/2.0_rt;
347 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
349 const amrex::Real s_x = (xint < 0.5_rt) ?
350 -1.0_rt/2.0_rt*xint_3 - 1.0_rt/4.0_rt*xint_2 + (5.0_rt/8.0_rt)*xint + 23.0_rt/48.0_rt :
351 (1.0_rt/2.0_rt)*xint_3 - 7.0_rt/4.0_rt*xint_2 + (11.0_rt/8.0_rt)*xint + 17.0_rt/48.0_rt;
352 const amrex::Real sdx = xint_3 - 3.0_rt/2.0_rt*xint_2 - 1.0_rt/2.0_rt*xint + 1.0_rt/2.0_rt;
353 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
355 const amrex::Real s_x = (xint < 0.5_rt) ?
356 (1.0_rt/6.0_rt)*xint_3 + (1.0_rt/4.0_rt)*xint_2 + (1.0_rt/8.0_rt)*xint + 1.0_rt/48.0_rt :
357 -1.0_rt/2.0_rt*xint_3 + (5.0_rt/4.0_rt)*xint_2 - 3.0_rt/8.0_rt*xint + 5.0_rt/48.0_rt;
358 const amrex::Real sdx = -2.0_rt/3.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 + (1.0_rt/2.0_rt)*xint + 1.0_rt/6.0_rt;
359 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
361 const amrex::Real s_x = (xint < 0.5_rt) ?
363 (1.0_rt/6.0_rt)*xint_3 - 1.0_rt/4.0_rt*xint_2 + (1.0_rt/8.0_rt)*xint - 1.0_rt/48.0_rt;
364 const amrex::Real sdx = (1.0_rt/6.0_rt)*xint_3;
365 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
368 }
else if constexpr (derivative_type==2) {
369 if constexpr (depos_order==0) {
371 const amrex::Real xfloor = amrex::Math::floor(xmid);
373 const amrex::Real s_x = 0;
374 const amrex::Real sdx = -1.0_rt/2.0_rt;
375 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
377 const amrex::Real s_x = 1;
378 const amrex::Real sdx = 0;
379 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
381 const amrex::Real s_x = 0;
382 const amrex::Real sdx = 1.0_rt/2.0_rt;
383 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
385 }
else if constexpr (depos_order==1) {
386 const amrex::Real xfloor = amrex::Math::floor(xmid);
387 const amrex::Real xint = xmid-xfloor;
389 const amrex::Real s_x = 0;
390 const amrex::Real sdx = (1.0_rt/2.0_rt)*xint - 1.0_rt/2.0_rt;
391 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
393 const amrex::Real s_x = 1 - xint;
394 const amrex::Real sdx = -1.0_rt/2.0_rt*xint;
395 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
397 const amrex::Real s_x = xint;
398 const amrex::Real sdx = 1.0_rt/2.0_rt - 1.0_rt/2.0_rt*xint;
399 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
401 const amrex::Real s_x = 0;
402 const amrex::Real sdx = (1.0_rt/2.0_rt)*xint;
403 return {s_x, -sdx,
static_cast<int>(xfloor)-1+ix};
405 }
else if constexpr (depos_order==2) {
407 const amrex::Real xfloor = amrex::Math::floor(xmid);
408 const amrex::Real xint = xmid-xfloor;
409 const amrex::Real xint_2 = xint*xint;
411 const amrex::Real s_x = 0;
412 const amrex::Real sdx = -1.0_rt/4.0_rt*xint_2 + (1.0_rt/2.0_rt)*xint - 1.0_rt/4.0_rt;
413 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
415 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_2 - xint + 1.0_rt/2.0_rt;
416 const amrex::Real sdx = (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint - 1.0_rt/4.0_rt;
417 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
419 const amrex::Real s_x = -xint_2 + xint + 1.0_rt/2.0_rt;
420 const amrex::Real sdx = 1.0_rt/4.0_rt - 1.0_rt/2.0_rt*xint;
421 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
423 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_2;
424 const amrex::Real sdx = -1.0_rt/2.0_rt*xint_2 + (1.0_rt/2.0_rt)*xint + 1.0_rt/4.0_rt;
425 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
427 const amrex::Real s_x = 0;
428 const amrex::Real sdx = (1.0_rt/4.0_rt)*xint_2;
429 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
431 }
else if constexpr (depos_order==3) {
432 const amrex::Real xfloor = amrex::Math::floor(xmid);
433 const amrex::Real xint = xmid-xfloor;
434 const amrex::Real xint_2 = xint*xint;
435 const amrex::Real xint_3 = xint_2*xint;
437 const amrex::Real s_x = 0;
438 const amrex::Real sdx = (1.0_rt/12.0_rt)*xint_3 - 1.0_rt/4.0_rt*xint_2 + (1.0_rt/4.0_rt)*xint - 1.0_rt/12.0_rt;
439 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
441 const amrex::Real s_x = -1.0_rt/6.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/2.0_rt*xint + 1.0_rt/6.0_rt;
442 const amrex::Real sdx = -1.0_rt/4.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 - 1.0_rt/3.0_rt;
443 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
445 const amrex::Real s_x = (1.0_rt/2.0_rt)*xint_3 - xint_2 + 2.0_rt/3.0_rt;
446 const amrex::Real sdx = (1.0_rt/6.0_rt)*xint_3 - 1.0_rt/2.0_rt*xint;
447 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
449 const amrex::Real s_x = -1.0_rt/2.0_rt*xint_3 + (1.0_rt/2.0_rt)*xint_2 + (1.0_rt/2.0_rt)*xint + 1.0_rt/6.0_rt;
450 const amrex::Real sdx = (1.0_rt/6.0_rt)*xint_3 - 1.0_rt/2.0_rt*xint_2 + 1.0_rt/3.0_rt;
451 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
453 const amrex::Real s_x = (1.0_rt/6.0_rt)*xint_3;
454 const amrex::Real sdx = -1.0_rt/4.0_rt*xint_3 + (1.0_rt/4.0_rt)*xint_2 + (1.0_rt/4.0_rt)*xint + 1.0_rt/12.0_rt;
455 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
457 const amrex::Real s_x = 0;
458 const amrex::Real sdx = (1.0_rt/12.0_rt)*xint_3;
459 return {s_x, -sdx,
static_cast<int>(xfloor)-2+ix};
463 static_assert(0 <= depos_order && depos_order <= 3 &&
464 0 <= derivative_type && derivative_type <= 2);
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int compute_shape_factor< 0 >(amrex::Real *const sx, amrex::Real xmid)
Definition: ShapeFactors.H:42
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE derivative_shape_factor_result single_derivative_shape_factor(amrex::Real xmid, int ix) noexcept
Compute a single derivative shape factor and return the index of the cell where the particle writes.
Definition: ShapeFactors.H:212
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE shape_factor_result compute_single_shape_factor(amrex::Real xmid, int ix) noexcept
Compute a single shape factor and return the index of the cell where the particle writes.
Definition: ShapeFactors.H:123
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int compute_shape_factor< 2 >(amrex::Real *const sx, amrex::Real xmid)
Definition: ShapeFactors.H:76
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int compute_shape_factor< 1 >(amrex::Real *const sx, amrex::Real xmid)
Definition: ShapeFactors.H:58
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int compute_shape_factor< 3 >(amrex::Real *const sx, amrex::Real xmid)
Definition: ShapeFactors.H:96
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int compute_shape_factor(amrex::Real *const sx, amrex::Real xint)
Definition: ShapeFactors.H:29
j
Definition: MakeOpenBoundary.py:128
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
Definition: ShapeFactors.H:197
int cell
Definition: ShapeFactors.H:200
amrex::Real factor
Definition: ShapeFactors.H:198
amrex::Real dx_factor
Definition: ShapeFactors.H:199
Definition: ShapeFactors.H:110
amrex::Real factor
Definition: ShapeFactors.H:111
int cell
Definition: ShapeFactors.H:112