GInX
Geodesics Integrator tool for particles in GRMHD using Adaptive Mesh Refinement using AMReX.
Interpolator.hxx
Go to the documentation of this file.
1
8#ifndef INTERPOLATOR_HXX
9#define INTERPOLATOR_HXX
10
11#include "AMReX_Box.H"
12#include "AMReX_Config.H"
13#include "AMReX_MFIter.H"
14#include "AMReX_REAL.H"
15#include "AMReX_Random.H"
16#include "AMReX_RandomEngine.H"
17#include "AMReX_Scan.H"
18#include "cctk_Types.h"
19#include <array>
20#include <cctk.h>
21#include <iostream>
22
26namespace GInX {
27
28// #############################################################################
29// Barycentric Lagrange Interpolator
30// #############################################################################
31
52template <int N>
53AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
54barycentric_cubic_1d(CCTK_REAL &value, const int (&points)[N],
55 const CCTK_REAL *weights, const CCTK_REAL (&values)[N],
56 const CCTK_REAL &x, const CCTK_REAL &plo,
57 const CCTK_REAL &dx) {
58 CCTK_REAL num{0.0};
59 CCTK_REAL den{0.0};
60
61 for (int i = 0; i < N; i++) {
62 // Check if the point belongs to the consulted points.
63 CCTK_REAL diff = x - (plo + points[i] * dx);
64 const CCTK_REAL tolerance = 1e-12;
65 if (diff < tolerance && diff > -tolerance) {
66 value = values[i];
67 return;
68 }
69 // Compute the weights and values
70 CCTK_REAL term = weights[i] / diff;
71 num += term * values[i];
72 den += term;
73 }
74
75 // Return the interpolation
76 value = num / den;
77}
78
95template <int INTERPOLATION_ORDER>
96AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
97barycentric_cubic_3d(const amrex::Array4<CCTK_REAL const> &f,
98 const long int &i0, const long int &j0, const long int &k0,
99 const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z,
100 const amrex::GpuArray<double, 3> &dx,
101 const amrex::GpuArray<double, 3> &plo, const int &comp) {
102 const int order =
103 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
104
105 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
106
107 const CCTK_REAL all_weights[14] = {
108 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
109 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
110 const CCTK_INT *nodes = &all_nodes[order];
111 const CCTK_REAL *w = &all_weights[order];
112
113 // Do the interpolation on x
114 CCTK_REAL G[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
115 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
116 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
117 CCTK_REAL values[INTERPOLATION_ORDER];
118 CCTK_INT points[INTERPOLATION_ORDER];
119 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
120 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], comp);
121 points[i] = i0 + nodes[i];
122 }
123 barycentric_cubic_1d<INTERPOLATION_ORDER>(G[j][k], points, w, values, x,
124 plo[0], dx[0]);
125 }
126 }
127
128 // Do the interpolation on y
129 CCTK_REAL H[INTERPOLATION_ORDER];
130 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
131 CCTK_REAL values[INTERPOLATION_ORDER];
132 CCTK_INT points[INTERPOLATION_ORDER];
133 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
134 values[j] = G[j][k];
135 points[j] = j0 + nodes[j];
136 }
137 barycentric_cubic_1d<INTERPOLATION_ORDER>(H[k], points, w, values, y,
138 plo[1], dx[1]);
139 }
140
141 // Do the interpolation on z
142 CCTK_INT points[INTERPOLATION_ORDER];
143 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
144 points[k] = k0 + nodes[k];
145 }
146 CCTK_REAL value;
147 barycentric_cubic_1d<INTERPOLATION_ORDER>(value, points, w, H, z, plo[2],
148 dx[2]);
149 return value;
150} // barycentric_cubic_3d with component
151
187template <int INTERPOLATION_ORDER>
188AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL
189barycentric_cubic_3d(const amrex::Array4<CCTK_REAL const> &f,
190 const long int &i0, const long int &j0, const long int &k0,
191 const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z,
192 const amrex::GpuArray<double, 3> &dx,
193 const amrex::GpuArray<double, 3> &plo) {
194
195 return barycentric_cubic_3d<INTERPOLATION_ORDER>(f, i0, j0, k0, x, y, z, dx,
196 plo, 0);
197} // barycentric_cubic_3d No component
198
242template <int N>
243AMREX_GPU_DEVICE
244 AMREX_GPU_HOST AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
245 der_barycentric_cubic_1d(CCTK_REAL &f_x, CCTK_REAL &d_f_x,
246 const int (&points)[N], const CCTK_REAL *weights,
247 const CCTK_REAL (&values)[N], const CCTK_REAL &x,
248 const CCTK_REAL &plo, const CCTK_REAL &dx) {
249 CCTK_REAL num{0.0};
250 CCTK_REAL den{0.0};
251 CCTK_REAL den_sqr{0.0};
252 CCTK_REAL der_num{0.0};
253 d_f_x = 0.0;
254
255 // Compute the interpolation
256 for (int i = 0; i < N; i++) {
257 CCTK_REAL diff = x - (plo + points[i] * dx);
258 const CCTK_REAL tolerance = 1e-12;
259 if (diff < tolerance && diff > -tolerance) {
260 // Check if the point makes part of the points used on the
261 // interpolation.
262 for (int j = 0; j < N; j++) {
263 if (i == j) {
264 continue;
265 }
266 const auto x_j = points[j] * dx;
267 const auto x_i = points[i] * dx;
268 d_f_x += weights[j] * (values[i] - values[j]) / (x_j - x_i);
269 }
270 d_f_x /= weights[i];
271 f_x = values[i];
272 return;
273 }
274
275 // Compute the weights for the interpolation and the derivative.
276 CCTK_REAL term = weights[i] / diff;
277 num += term * values[i];
278 den += term;
279 den_sqr += term / diff;
280 }
281
282 // Compute the weights for the derivative.
283 for (int i = 0; i < N; i++) {
284 CCTK_REAL term = weights[i] / (x - (plo + points[i] * dx));
285 der_num += (-term * den / (x - (plo + points[i] * dx)) + term * den_sqr) *
286 values[i];
287 }
288
289 // Fill the values
290 f_x = num / den;
291 d_f_x = der_num / (den * den);
292} // der_barycentric_cubic_1d
293
363template <int INTERPOLATION_ORDER>
364AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
365barycentric_derivative_and_interpolate(CCTK_REAL &f_xyz, CCTK_REAL &df_xyz_0,
366 CCTK_REAL &df_xyz_1, CCTK_REAL &df_xyz_2,
367 amrex::Array4<CCTK_REAL const> const &f,
368 const long int &i0, const long int &j0,
369 const long int &k0, const CCTK_REAL &x,
370 const CCTK_REAL &y, const CCTK_REAL &z,
371 const amrex::GpuArray<double, 3> &dx,
372 const amrex::GpuArray<double, 3> &plo,
373 const int &comp) {
374
375 const int order =
376 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
377 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
378 const CCTK_REAL all_weights[14] = {
379 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
380 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
381 const CCTK_INT *nodes = &all_nodes[order];
382 const CCTK_REAL *w = &all_weights[order];
383
384 // Computing f(x, y_i, z_i)
385 CCTK_REAL G_xyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
386 // Computing d/dx f(x, y_i, z_i)
387 CCTK_REAL G_dxyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
388 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
389 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
390 CCTK_REAL values[INTERPOLATION_ORDER];
391 CCTK_INT points[INTERPOLATION_ORDER];
392 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
393 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], comp);
394 points[i] = i0 + nodes[i];
395 }
396 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
397 G_xyz[j][k], G_dxyz[j][k], points, w, values, x, plo[0], dx[0]);
398 }
399 }
400
401 // Computing f(x, y, z_i)
402 CCTK_REAL H_xyz[INTERPOLATION_ORDER];
403 // Computing d/dx f(x, y, z_i)
404 CCTK_REAL H_dxyz[INTERPOLATION_ORDER];
405 // Computing d/dy f(x, y, z_i)
406 CCTK_REAL H_xdyz[INTERPOLATION_ORDER];
407 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
408 CCTK_REAL values[INTERPOLATION_ORDER];
409 CCTK_REAL d_values[INTERPOLATION_ORDER];
410 CCTK_INT points[INTERPOLATION_ORDER];
411 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
412 values[j] = G_xyz[j][k];
413 d_values[j] = G_dxyz[j][k];
414 points[j] = j0 + nodes[j];
415 }
416 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(H_xyz[k], H_xdyz[k], points,
417 w, values, y, plo[1], dx[1]);
418 barycentric_cubic_1d<INTERPOLATION_ORDER>(H_dxyz[k], points, w, d_values, y,
419 plo[1], dx[1]);
420 }
421
422 CCTK_INT points[INTERPOLATION_ORDER];
423 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
424 points[k] = k0 + nodes[k];
425 }
426 // Computing f(x, y, z)
427 // Computing d/dz f(x, y, z)
428 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(f_xyz, df_xyz_2, points, w,
429 H_xyz, z, plo[2], dx[2]);
430 // Computing d/dx f(x, y, z)
431 barycentric_cubic_1d<INTERPOLATION_ORDER>(df_xyz_0, points, w, H_dxyz, z,
432 plo[2], dx[2]);
433 // Computing d/dy f(x, y, z)
434 barycentric_cubic_1d<INTERPOLATION_ORDER>(df_xyz_1, points, w, H_xdyz, z,
435 plo[2], dx[2]);
436} // barycentric_derivative_and_interpolate
437
462template <int INTERPOLATION_ORDER>
463AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
464interpolate_array(amrex::GpuArray<CCTK_REAL, 6> &array6,
465 const amrex::Array4<CCTK_REAL const> &f, const long int &i0,
466 const long int &j0, const long int &k0, const CCTK_REAL &x,
467 const CCTK_REAL &y, const CCTK_REAL &z,
468 const amrex::GpuArray<double, 3> &dx,
469 const amrex::GpuArray<double, 3> &plo) {
470 const int order =
471 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
472 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
473 const CCTK_REAL all_weights[14] = {
474 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
475 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
476 const CCTK_INT *nodes = &all_nodes[order];
477 const CCTK_REAL *w = &all_weights[order];
478
479 for (int comp = 0; comp < 6; comp++) {
480
481 // Do the interpolation on x
482 CCTK_REAL G[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
483 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
484 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
485 CCTK_REAL values[INTERPOLATION_ORDER];
486 CCTK_INT points[INTERPOLATION_ORDER];
487 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
488 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], comp);
489 points[i] = i0 + nodes[i];
490 }
491 barycentric_cubic_1d<INTERPOLATION_ORDER>(G[j][k], points, w, values, x,
492 plo[0], dx[0]);
493 }
494 }
495
496 // Do the interpolation on y
497 CCTK_REAL H[INTERPOLATION_ORDER];
498 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
499 CCTK_REAL values[INTERPOLATION_ORDER];
500 CCTK_INT points[INTERPOLATION_ORDER];
501 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
502 values[j] = G[j][k];
503 points[j] = j0 + nodes[j];
504 }
505 barycentric_cubic_1d<INTERPOLATION_ORDER>(H[k], points, w, values, y,
506 plo[1], dx[1]);
507 }
508
509 // Do the interpolation on z
510 CCTK_INT points[INTERPOLATION_ORDER];
511 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
512 points[k] = k0 + nodes[k];
513 }
514 barycentric_cubic_1d<INTERPOLATION_ORDER>(array6[comp], points, w, H, z,
515 plo[2], dx[2]);
516 }
517} // interpolate_array
518
545template <int INTERPOLATION_ORDER>
546AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
547d_interpolate_array(amrex::GpuArray<CCTK_REAL, 6> &array6,
548 amrex::GpuArray<amrex::GpuArray<CCTK_REAL, 6>, 3> &d_array6,
549 const amrex::Array4<CCTK_REAL const> &f, const long int &i0,
550 const long int &j0, const long int &k0, const CCTK_REAL &x,
551 const CCTK_REAL &y, const CCTK_REAL &z,
552 const amrex::GpuArray<double, 3> &dx,
553 const amrex::GpuArray<double, 3> &plo) {
554 const int order =
555 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
556 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
557 const CCTK_REAL all_weights[14] = {
558 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
559 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
560 const CCTK_INT *nodes = &all_nodes[order];
561 const CCTK_REAL *w = &all_weights[order];
562
563 for (int comp = 0; comp < 6; comp++) {
564
565 // Computing f(x, y_i, z_i)
566 CCTK_REAL G_xyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
567 // Computing d/dx f(x, y_i, z_i)
568 CCTK_REAL G_dxyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
569 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
570 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
571 CCTK_REAL values[INTERPOLATION_ORDER];
572 CCTK_INT points[INTERPOLATION_ORDER];
573 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
574 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], comp);
575 points[i] = i0 + nodes[i];
576 }
577 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
578 G_xyz[j][k], G_dxyz[j][k], points, w, values, x, plo[0], dx[0]);
579 }
580 }
581
582 // Computing f(x, y, z_i)
583 CCTK_REAL H_xyz[INTERPOLATION_ORDER];
584 // Computing d/dx f(x, y, z_i)
585 CCTK_REAL H_dxyz[INTERPOLATION_ORDER];
586 // Computing d/dy f(x, y, z_i)
587 CCTK_REAL H_xdyz[INTERPOLATION_ORDER];
588 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
589 CCTK_REAL values[INTERPOLATION_ORDER];
590 CCTK_REAL d_values[INTERPOLATION_ORDER];
591 CCTK_INT points[INTERPOLATION_ORDER];
592 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
593 values[j] = G_xyz[j][k];
594 d_values[j] = G_dxyz[j][k];
595 points[j] = j0 + nodes[j];
596 }
597 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
598 H_xyz[k], H_xdyz[k], points, w, values, y, plo[1], dx[1]);
599 barycentric_cubic_1d<INTERPOLATION_ORDER>(H_dxyz[k], points, w, d_values,
600 y, plo[1], dx[1]);
601 }
602
603 CCTK_INT points[INTERPOLATION_ORDER];
604 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
605 points[k] = k0 + nodes[k];
606 }
607 // Computing f(x, y, z)
608 // Computing d/dz f(x, y, z)
609 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
610 array6[comp], d_array6[2][comp], points, w, H_xyz, z, plo[2], dx[2]);
611 // Computing d/dx f(x, y, z)
612 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array6[0][comp], points, w,
613 H_dxyz, z, plo[2], dx[2]);
614 // Computing d/dy f(x, y, z)
615 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array6[1][comp], points, w,
616 H_xdyz, z, plo[2], dx[2]);
617 }
618} // d_interpolate_array
619
646template <int INTERPOLATION_ORDER>
647AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
648d_interpolate_array(amrex::GpuArray<CCTK_REAL, 3> &array3,
649 amrex::GpuArray<amrex::GpuArray<CCTK_REAL, 3>, 3> &d_array3,
650 const amrex::Array4<CCTK_REAL const> &f, const long int &i0,
651 const long int &j0, const long int &k0, const CCTK_REAL &x,
652 const CCTK_REAL &y, const CCTK_REAL &z,
653 const amrex::GpuArray<double, 3> &dx,
654 const amrex::GpuArray<double, 3> &plo) {
655 const int order =
656 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
657 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
658 const CCTK_REAL all_weights[14] = {
659 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
660 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
661 const CCTK_INT *nodes = &all_nodes[order];
662 const CCTK_REAL *w = &all_weights[order];
663
664 for (int comp = 0; comp < 3; comp++) {
665
666 // Computing f(x, y_i, z_i)
667 CCTK_REAL G_xyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
668 // Computing d/dx f(x, y_i, z_i)
669 CCTK_REAL G_dxyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
670 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
671 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
672 CCTK_REAL values[INTERPOLATION_ORDER];
673 CCTK_INT points[INTERPOLATION_ORDER];
674 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
675 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], comp);
676 points[i] = i0 + nodes[i];
677 }
678 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
679 G_xyz[j][k], G_dxyz[j][k], points, w, values, x, plo[0], dx[0]);
680 }
681 }
682
683 // Computing f(x, y, z_i)
684 CCTK_REAL H_xyz[INTERPOLATION_ORDER];
685 // Computing d/dx f(x, y, z_i)
686 CCTK_REAL H_dxyz[INTERPOLATION_ORDER];
687 // Computing d/dy f(x, y, z_i)
688 CCTK_REAL H_xdyz[INTERPOLATION_ORDER];
689 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
690 CCTK_REAL values[INTERPOLATION_ORDER];
691 CCTK_REAL d_values[INTERPOLATION_ORDER];
692 CCTK_INT points[INTERPOLATION_ORDER];
693 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
694 values[j] = G_xyz[j][k];
695 d_values[j] = G_dxyz[j][k];
696 points[j] = j0 + nodes[j];
697 }
698 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
699 H_xyz[k], H_xdyz[k], points, w, values, y, plo[1], dx[1]);
700 barycentric_cubic_1d<INTERPOLATION_ORDER>(H_dxyz[k], points, w, d_values,
701 y, plo[1], dx[1]);
702 }
703
704 CCTK_INT points[INTERPOLATION_ORDER];
705 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
706 points[k] = k0 + nodes[k];
707 }
708 // Computing f(x, y, z)
709 // Computing d/dz f(x, y, z)
710 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
711 array3[comp], d_array3[2][comp], points, w, H_xyz, z, plo[2], dx[2]);
712 // Computing d/dx f(x, y, z)
713 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array3[0][comp], points, w,
714 H_dxyz, z, plo[2], dx[2]);
715 // Computing d/dy f(x, y, z)
716 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array3[1][comp], points, w,
717 H_xdyz, z, plo[2], dx[2]);
718 }
719} // d_interpolate_array with size 3 arrays
720
747template <int INTERPOLATION_ORDER>
748AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void
749d_interpolate_array(CCTK_REAL &array, amrex::GpuArray<CCTK_REAL, 3> &d_array,
750 const amrex::Array4<CCTK_REAL const> &f, const long int &i0,
751 const long int &j0, const long int &k0, const CCTK_REAL &x,
752 const CCTK_REAL &y, const CCTK_REAL &z,
753 const amrex::GpuArray<double, 3> &dx,
754 const amrex::GpuArray<double, 3> &plo) {
755 const int order =
756 (((INTERPOLATION_ORDER - 1) * INTERPOLATION_ORDER) >> 1) - 1;
757 const CCTK_INT all_nodes[14] = {0, 1, -1, 0, 1, -1, 0, 1, 2, -2, -1, 0, 1, 2};
758 const CCTK_REAL all_weights[14] = {
759 -1., 1., 0.5, -1., 0.5, -1.0 / 6.0, 0.5,
760 -0.5, 1.0 / 6.0, 1.0 / 24.0, -1.0 / 6.0, 0.25, -1.0 / 6.0, 1.0 / 24.0};
761 const CCTK_INT *nodes = &all_nodes[order];
762 const CCTK_REAL *w = &all_weights[order];
763
764 // Computing f(x, y_i, z_i)
765 CCTK_REAL G_xyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
766 // Computing d/dx f(x, y_i, z_i)
767 CCTK_REAL G_dxyz[INTERPOLATION_ORDER][INTERPOLATION_ORDER];
768 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
769 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
770 CCTK_REAL values[INTERPOLATION_ORDER];
771 CCTK_INT points[INTERPOLATION_ORDER];
772 for (int i = 0; i < INTERPOLATION_ORDER; i++) {
773 values[i] = f(i0 + nodes[i], j0 + nodes[j], k0 + nodes[k], 0);
774 points[i] = i0 + nodes[i];
775 }
776 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(
777 G_xyz[j][k], G_dxyz[j][k], points, w, values, x, plo[0], dx[0]);
778 }
779 }
780
781 // Computing f(x, y, z_i)
782 CCTK_REAL H_xyz[INTERPOLATION_ORDER];
783 // Computing d/dx f(x, y, z_i)
784 CCTK_REAL H_dxyz[INTERPOLATION_ORDER];
785 // Computing d/dy f(x, y, z_i)
786 CCTK_REAL H_xdyz[INTERPOLATION_ORDER];
787 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
788 CCTK_REAL values[INTERPOLATION_ORDER];
789 CCTK_REAL d_values[INTERPOLATION_ORDER];
790 CCTK_INT points[INTERPOLATION_ORDER];
791 for (int j = 0; j < INTERPOLATION_ORDER; j++) {
792 values[j] = G_xyz[j][k];
793 d_values[j] = G_dxyz[j][k];
794 points[j] = j0 + nodes[j];
795 }
796 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(H_xyz[k], H_xdyz[k], points,
797 w, values, y, plo[1], dx[1]);
798 barycentric_cubic_1d<INTERPOLATION_ORDER>(H_dxyz[k], points, w, d_values, y,
799 plo[1], dx[1]);
800 }
801
802 CCTK_INT points[INTERPOLATION_ORDER];
803 for (int k = 0; k < INTERPOLATION_ORDER; k++) {
804 points[k] = k0 + nodes[k];
805 }
806 // Computing f(x, y, z)
807 // Computing d/dz f(x, y, z)
808 der_barycentric_cubic_1d<INTERPOLATION_ORDER>(array, d_array[2], points, w,
809 H_xyz, z, plo[2], dx[2]);
810 // Computing d/dx f(x, y, z)
811 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array[0], points, w, H_dxyz, z,
812 plo[2], dx[2]);
813 // Computing d/dy f(x, y, z)
814 barycentric_cubic_1d<INTERPOLATION_ORDER>(d_array[1], points, w, H_xdyz, z,
815 plo[2], dx[2]);
816} // d_interpolate_array scalar
817
818} // namespace GInX
819
820#endif // !INTERPOLATOR_HXX
GInX
Interpolators namespace.
Definition: ParticlesContainer.hxx:39
GInX::der_barycentric_cubic_1d
AMREX_GPU_DEVICE AMREX_GPU_HOST AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void der_barycentric_cubic_1d(CCTK_REAL &f_x, CCTK_REAL &d_f_x, const int(&points)[N], const CCTK_REAL *weights, const CCTK_REAL(&values)[N], const CCTK_REAL &x, const CCTK_REAL &plo, const CCTK_REAL &dx)
Do a Barycentric interpolation and its derivative in one direction.
Definition: Interpolator.hxx:245
GInX::interpolate_array
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void interpolate_array(amrex::GpuArray< CCTK_REAL, 6 > &array6, const amrex::Array4< CCTK_REAL const > &f, const long int &i0, const long int &j0, const long int &k0, const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z, const amrex::GpuArray< double, 3 > &dx, const amrex::GpuArray< double, 3 > &plo)
Do the barycentric interpolation for a 3 dimensional and symmetric matrix.
Definition: Interpolator.hxx:464
GInX::d_interpolate_array
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void d_interpolate_array(amrex::GpuArray< CCTK_REAL, 6 > &array6, amrex::GpuArray< amrex::GpuArray< CCTK_REAL, 6 >, 3 > &d_array6, const amrex::Array4< CCTK_REAL const > &f, const long int &i0, const long int &j0, const long int &k0, const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z, const amrex::GpuArray< double, 3 > &dx, const amrex::GpuArray< double, 3 > &plo)
Do the barycentric interpolation for a 3 dimensional and symmetric matrix and its derivatives.
Definition: Interpolator.hxx:547
GInX::barycentric_cubic_1d
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void barycentric_cubic_1d(CCTK_REAL &value, const int(&points)[N], const CCTK_REAL *weights, const CCTK_REAL(&values)[N], const CCTK_REAL &x, const CCTK_REAL &plo, const CCTK_REAL &dx)
Do a Barycentric interpolation in one direction.
Definition: Interpolator.hxx:54
GInX::barycentric_cubic_3d
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_REAL barycentric_cubic_3d(const amrex::Array4< CCTK_REAL const > &f, const long int &i0, const long int &j0, const long int &k0, const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z, const amrex::GpuArray< double, 3 > &dx, const amrex::GpuArray< double, 3 > &plo, const int &comp)
Do a Barycentric interpolation in three direction for a vectorial function.
Definition: Interpolator.hxx:97
GInX::barycentric_derivative_and_interpolate
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CCTK_ATTRIBUTE_ALWAYS_INLINE void barycentric_derivative_and_interpolate(CCTK_REAL &f_xyz, CCTK_REAL &df_xyz_0, CCTK_REAL &df_xyz_1, CCTK_REAL &df_xyz_2, amrex::Array4< CCTK_REAL const > const &f, const long int &i0, const long int &j0, const long int &k0, const CCTK_REAL &x, const CCTK_REAL &y, const CCTK_REAL &z, const amrex::GpuArray< double, 3 > &dx, const amrex::GpuArray< double, 3 > &plo, const int &comp)
Do a Barycentric interpolation and its derivatives in three directions.
Definition: Interpolator.hxx:365