eschnett/cactusamrex - Change IC2NW4EHRU42VTEF6RNOC26DVB3BKSA6JLGMO6CUXGIHHOZ3KTCQC

CarpetX: Support higher interpolation orders

Created by Erik Schnetter on February 27, 2020

IC2NW4EHRU42VTEF6RNOC26DVB3BKSA6JLGMO6CUXGIHHOZ3KTCQC

Dependencies

In channels

master

main

Change contents

Insertion in CarpetX/param.ccl at line 201 [5.17918]

[3.1]

[4.0]

CCTK_INT interpolation_order "Interpolation order" STEERABLE=always
{
  0:* :: ""
} 1

Insertion in CarpetX/param.ccl at line 207 [5.17918]
[3.72]
[4.1]

Deletion in CarpetX/src/interpolate.cxx at line 25 [2.373]

B:BD[2.788] → [2.788:914]

// Calculate 10^n
constexpr inline int iexp10(int n) {
  int r = 1;
  while (n > 0) {
    r *= 10;
    --n;
  }
  return r;
}

Insertion in CarpetX/src/interpolate.cxx at line 48 [2.373]

[2.1662]

    case 0: {
      const T x = di[dir];
#ifdef CCTK_DEBUG
      assert(fabs(x) <= T(0.5));
#endif
      const T y0 = interpolate<dir - 1>(i, di);
      switch (derivs[dir]) {
      case 0:
        return y0;
      case 1:
        return 0;
      case 2:
        return 0;
      }
    }

Replacement in CarpetX/src/interpolate.cxx at line 64 [2.373]

B:BD[2.1676] → [2.1676:1765]

      T y0 = interpolate<dir - 1>(i, di);
      T y1 = interpolate<dir - 1>(i + DI, di);

[2.1676]

[2.1765]

      const T x = di[dir];
#ifdef CCTK_DEBUG
      assert(x >= T(0) && x <= T(1));
#endif
      const T y0 = interpolate<dir - 1>(i, di);
      const T y1 = interpolate<dir - 1>(i + DI, di);

Replacement in CarpetX/src/interpolate.cxx at line 72 [2.373]
B:BD[2.1808] → [2.1808:1858]
```
        return (1 - di[dir]) * y0 + di[dir] * y1;
```
[2.1808]
[2.1858]
```
        return (1 - x) * y0 + x * y1;
```
Insertion in CarpetX/src/interpolate.cxx at line 75 [2.373]
[2.1909]
[2.1909]
```
      case 2:
        return 0;
```

Replacement in CarpetX/src/interpolate.cxx at line 80 [2.373]

B:BD[2.1937] → [2.1937:2077]

      T y0 = interpolate<dir - 1>(i, di);
      T y1 = interpolate<dir - 1>(i + DI, di);
      T y2 = interpolate<dir - 1>(i + 2 * DI, di);

[2.1937]

[2.2077]

      const T x = di[dir] - order / T(2);
#ifdef CCTK_DEBUG
      assert(fabs(x) <= T(0.5));
#endif
      const T y0 = interpolate<dir - 1>(i, di);
      const T y1 = interpolate<dir - 1>(i + DI, di);
      const T y2 = interpolate<dir - 1>(i + 2 * DI, di);

Replacement in CarpetX/src/interpolate.cxx at line 89 [2.373]

B:BD[2.2120] → [2.2120:2276]

        return (di[dir] - 1) * (di[dir] - 2) / 2 * y0 +
               di[dir] * (di[dir] - 2) / -1 * y1 +
               di[dir] * (di[dir] - 1) / 2 * y2;

[2.2120]

[2.2276]

        return (-1 / T(2) * x + 1 / T(2) * pow(x, 2)) * y0 +
               (1 - pow(x, 2)) * y1 +
               (1 / T(2) * x + 1 / T(2) * pow(x, 2)) * y2;

Replacement in CarpetX/src/interpolate.cxx at line 93 [2.373]

B:BD[2.2290] → [2.2290:2411]

        return ((2 * di[dir] - 3) / 2 * y0 + (2 * di[dir] - 2) / -1 * y1 +
                (2 * di[dir] - 1) / 2 * y2) /

[2.2290]

[2.2411]

        return ((-1 / T(2) + x) * y0 - 2 * x * y1 + (1 / T(2) + x) * y2) /

Insertion in CarpetX/src/interpolate.cxx at line 97 [2.373]

[2.2502]

      }
    }
    case 3: {
      const T x = di[dir] - order / T(2);
#ifdef CCTK_DEBUG
      assert(fabs(x) <= T(0.5));
#endif
      const T y0 = interpolate<dir - 1>(i, di);
      const T y1 = interpolate<dir - 1>(i + DI, di);
      const T y2 = interpolate<dir - 1>(i + 2 * DI, di);
      const T y3 = interpolate<dir - 1>(i + 3 * DI, di);
      switch (derivs[dir]) {
      case 0:
        return (-1 / T(16) + 1 / T(24) * x + 1 / T(4) * pow(x, 2) -
                1 / T(6) * pow(x, 3)) *
                   y0 +
               (9 / T(16) - 9 / T(8) * x - 1 / T(4) * pow(x, 2) +
                1 / T(2) * pow(x, 3)) *
                   y1 +
               (9 / T(16) + 9 / T(8) * x - 1 / T(4) * pow(x, 2) -
                1 / T(2) * pow(x, 3)) *
                   y2 +
               (-1 / T(16) - 1 / T(24) * x + 1 / T(4) * pow(x, 2) +
                1 / T(6) * pow(x, 3)) *
                   y3;
      case 1:
        return ((1 / T(24) + 1 / T(2) * x - 1 / T(2) * pow(x, 2)) * y0 +
                (-9 / T(8) - 1 / T(2) * x + 3 / T(2) * pow(x, 2)) * y1 +
                (9 / T(8) - 1 / T(2) * x - 3 / T(2) * pow(x, 2)) * y2 +
                (-1 / T(24) + 1 / T(2) * x + 1 / T(2) * pow(x, 2)) * y3) /
               dx[dir];
      case 2:
        return ((1 / T(2) - x) * y0 + (-1 / T(2) + 3 * x) * y1 +
                (-1 / T(2) - 3 * x) * y2 + (1 / T(2) + x) * y3) /
               pow(dx[dir], 2);

Insertion in CarpetX/src/interpolate.cxx at line 134 [2.373]

[2.2516]

    case 4: {
      const T x = di[dir] - order / T(2);
#ifdef CCTK_DEBUG
      assert(fabs(x) <= T(0.5));
#endif
      const T y0 = interpolate<dir - 1>(i, di);
      const T y1 = interpolate<dir - 1>(i + DI, di);
      const T y2 = interpolate<dir - 1>(i + 2 * DI, di);
      const T y3 = interpolate<dir - 1>(i + 3 * DI, di);
      const T y4 = interpolate<dir - 1>(i + 4 * DI, di);
      switch (derivs[dir]) {
      case 0:
        return (1 / T(12) * x - 1 / T(24) * pow(x, 2) - 1 / T(12) * pow(x, 3) +
                1 / T(24) * pow(x, 4)) *
                   y0 +
               (-2 / T(3) * x + 2 / T(3) * pow(x, 2) + 1 / T(6) * pow(x, 3) -
                1 / T(6) * pow(x, 4)) *
                   y1 +
               (1 - 5 / T(4) * pow(x, 2) + 1 / T(4) * pow(x, 4)) * y2 +
               (2 / T(3) * x + 2 / T(3) * pow(x, 2) - 1 / T(6) * pow(x, 3) -
                1 / T(6) * pow(x, 4)) *
                   y3 +
               (-1 / T(12) * x - 1 / T(24) * pow(x, 2) + 1 / T(12) * pow(x, 3) +
                1 / T(24) * pow(x, 4)) *
                   y4;
      case 1:
        return ((1 / T(12) - 1 / T(12) * x - 1 / T(4) * pow(x, 2) +
                 1 / T(6) * pow(x, 3)) *
                    y0 +
                (-2 / T(3) + 4 / T(3) * x + 1 / T(2) * pow(x, 2) -
                 2 / T(3) * pow(x, 3)) *
                    y1 +
                (-5 / T(2) * x + pow(x, 3)) * y2 +
                (2 / T(3) + 4 / T(3) * x - 1 / T(2) * pow(x, 2) -
                 2 / T(3) * pow(x, 3)) *
                    y3 +
                (-1 / T(12) - 1 / T(12) * x + 1 / T(4) * pow(x, 2) +
                 1 / T(6) * pow(x, 3)) *
                    y4) /
               dx[dir];
      case 2:
        return ((-1 / T(12) - 1 / T(2) * x + 1 / T(2) * pow(x, 2)) * y0 +
                (4 / T(3) + x - 2 * pow(x, 2)) * y1 +
                (-5 / T(2) + 3 * pow(x, 2)) * y2 +
                (4 / T(3) - x - 2 * pow(x, 2)) * y3 +
                (-1 / T(12) + 1 / T(2) * x + 1 / T(2) * pow(x, 2)) * y4) /
               pow(dx[dir], 2);
      }

Replacement in CarpetX/src/interpolate.cxx at line 183 [2.373]
B:BD[2.2522] → [2.2522:2537]
```
    assert(0);
```
[2.2522]
[2.2537]
```
    default:
      assert(0);
    }
```
Insertion in CarpetX/src/interpolate.cxx at line 200 [2.373]
[2.3208]
[2.3208]
```
  DECLARE_CCTK_PARAMETERS;
```
Deletion in CarpetX/src/interpolate.cxx at line 204 [2.373]
B:BD[2.3315] → [2.3315:3343]
```
  constexpr int order = 2;
```
Deletion in CarpetX/src/interpolate.cxx at line 273 [2.373]
B:BD[2.5730] → [2.5730:5768]
```
        const int op = operations[v];
```

Replacement in CarpetX/src/interpolate.cxx at line 274 [2.373]

B:BD[2.5799] → [2.5799:5882]

        for (int d = 0; d < dim; ++d)
          derivs[d] = (op / iexp10(d)) % 10;

[2.5799]

[2.5882]

        int op = operations[v];
        while (op > 0) {
          const int dir = op % 10 - 1;
          if (dir >= 0) {
            assert(dir >= 0 && dir < dim);
            ++derivs[dir];
          }
          op /= 10;
        }

Insertion in CarpetX/src/interpolate.cxx at line 286 [2.373]

[2.5957]

        switch (interpolation_order) {
        case 0: {
          constexpr int order = 0;
#pragma omp simd
          for (int n = 0; n < np; ++n) {
            vect<int, dim> i;
            vect<CCTK_REAL, dim> di;
            for (int d = 0; d < dim; ++d) {
              CCTK_REAL x = particles[n].pos(d);
              CCTK_REAL ri = (x - x0[d]) / dx[d];
              i[d] = lrint(ri - (order / CCTK_REAL(2)));
              di[d] = ri - i[d];
            }
            const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
            varresult.at(n) = interp.interpolate<dim - 1>(i, di);
          }
          break;
        }
        case 1: {
          constexpr int order = 1;

Replacement in CarpetX/src/interpolate.cxx at line 307 [2.373]

B:BD[2.5974] → [2.5974:6299]

        for (int n = 0; n < np; ++n) {
          vect<int, dim> i;
          vect<CCTK_REAL, dim> di;
          for (int d = 0; d < dim; ++d) {
            CCTK_REAL x = particles[n].pos(d);
            CCTK_REAL ri = (x - x0[d]) / dx[d];
            i[d] = lrint(ri - (order / CCTK_REAL(2)));
            di[d] = ri - i[d];

[2.5974]

[2.6299]

          for (int n = 0; n < np; ++n) {
            vect<int, dim> i;
            vect<CCTK_REAL, dim> di;
            for (int d = 0; d < dim; ++d) {
              CCTK_REAL x = particles[n].pos(d);
              CCTK_REAL ri = (x - x0[d]) / dx[d];
              i[d] = lrint(ri - (order / CCTK_REAL(2)));
              di[d] = ri - i[d];
            }
            const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
            varresult.at(n) = interp.interpolate<dim - 1>(i, di);
          }
          break;
        }
        case 2: {
          constexpr int order = 2;
#pragma omp simd
          for (int n = 0; n < np; ++n) {
            vect<int, dim> i;
            vect<CCTK_REAL, dim> di;
            for (int d = 0; d < dim; ++d) {
              CCTK_REAL x = particles[n].pos(d);
              CCTK_REAL ri = (x - x0[d]) / dx[d];
              i[d] = lrint(ri - (order / CCTK_REAL(2)));
              di[d] = ri - i[d];
            }
            const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
            varresult.at(n) = interp.interpolate<dim - 1>(i, di);
          }
          break;
        }
        case 3: {
          constexpr int order = 3;
#pragma omp simd
          for (int n = 0; n < np; ++n) {
            vect<int, dim> i;
            vect<CCTK_REAL, dim> di;
            for (int d = 0; d < dim; ++d) {
              CCTK_REAL x = particles[n].pos(d);
              CCTK_REAL ri = (x - x0[d]) / dx[d];
              i[d] = lrint(ri - (order / CCTK_REAL(2)));
              di[d] = ri - i[d];
            }
            const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
            varresult.at(n) = interp.interpolate<dim - 1>(i, di);

Replacement in CarpetX/src/interpolate.cxx at line 353 [2.373]

B:BD[2.6311] → [2.6311:7751]


          // const vect<int, dim> DI{1, 0, 0};
          // const vect<int, dim> DJ{0, 1, 0};
          // const vect<int, dim> DK{0, 0, 1};
          // const auto interp0{[&](const vect<int, dim> i) {
          //   return vars(i[0], i[1], i[2], vi);
          // }};
          // const auto interp1{[&](const vect<int, dim> i) {
          //   if (derivs[0] == 0)
          //     return (1 - di[0]) * interp0(i) + di[0] * interp0(i + DI);
          //   if (derivs[0] == 1)
          //     return (-interp0(i) + interp0(i + DI)) / dx[0];
          //   return 0 * interp0(i);
          // }};
          // const auto interp2{[&](const vect<int, dim> i) {
          //   if (derivs[1] == 0)
          //     return (1 - di[1]) * interp1(i) + di[1] * interp1(i + DJ);
          //   if (derivs[1] == 1)
          //     return (-interp1(i) + interp1(i + DJ)) / dx[1];
          //   return 0 * interp1(i);
          // }};
          // const auto interp3{[&](const vect<int, dim> i) {
          //   if (derivs[2] == 0)
          //     return (1 - di[2]) * interp2(i) + di[2] * interp2(i + DK);
          //   if (derivs[2] == 1)
          //     return (-interp2(i) + interp2(i + DK)) / dx[2];
          //   return 0 * interp2(i);
          // }};
          // varresult.at(n) = interp3(i);
          const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
          varresult.at(n) = interp.interpolate<dim - 1>(i, di);

[2.6311]

[2.7751]

          break;
        }
        case 4: {
          constexpr int order = 4;
#pragma omp simd
          for (int n = 0; n < np; ++n) {
            vect<int, dim> i;
            vect<CCTK_REAL, dim> di;
            for (int d = 0; d < dim; ++d) {
              CCTK_REAL x = particles[n].pos(d);
              CCTK_REAL ri = (x - x0[d]) / dx[d];
              i[d] = lrint(ri - (order / CCTK_REAL(2)));
              di[d] = ri - i[d];
            }
            const interpolator<CCTK_REAL, order> interp{vars, vi, derivs, dx};
            varresult.at(n) = interp.interpolate<dim - 1>(i, di);
          }
          break;
        }
        default:
          assert(0);