eschnett/cactusamrex - Change C7Y67X3DJBDM4QW3NAVSVJF6GSSLD466S4FFJSIKUEJANLPC6FIAC

CarpetX: Correct tests for conservative interpolation

Created by Erik Schnetter on January 4, 2020

C7Y67X3DJBDM4QW3NAVSVJF6GSSLD466S4FFJSIKUEJANLPC6FIAC

Dependencies

In channels

master

main

Change contents

Replacement in CarpetX/src/prolongate_3d_rf2.cxx at line 312 [4.1186]

B:BD[3.2495] → [3.2495:2688]

    for (int order = 0; order < ORDER; ++order) {
      auto f = [&](T x) { return order == 0 ? T(1) : pow(x, order); };
      auto fint = [&](T x) { return pow(x, order + 1) / (order + 1); };

[3.2495]

[3.2688]

    for (int order = 0; order <= ORDER; ++order) {
      // const auto f{[&](T x) {
      //   return (order + 1) * (order == 0 ? T(1) : pow(x, order));
      // }};
      const auto fint{[&](T x) { return pow(x, order + 1); }};

Replacement in CarpetX/src/prolongate_3d_rf2.cxx at line 323 [4.1186]

B:BD[3.2900] → [3.2900:2967]

          T x = (i - i0) + CENTERING / T(2);
          T y = f(x);

[3.2900]

[3.2967]

          const T x = (i - i0) + CENTERING / T(2);
          // T y = f(x);
          const T dx = 1;
          const T xlo = x - dx / 2;
          const T xhi = x + dx / 2;
          const T y = fint(xhi) - fint(xlo);

Replacement in CarpetX/src/prolongate_3d_rf2.cxx at line 339 [4.1186]

B:BD[3.3320] → [3.3320:3454]

        T dx = x1[1] - x1[0];
        T xlo = x1[0] - dx / 2;
        T xhi = x1[1] + dx / 2;
        T yint = fint(xhi) - fint(xlo);

[3.3320]

[3.3454]

        const T dx = x1[1] - x1[0];
        const T xlo = x1[0] - dx / 2;
        const T xhi = x1[1] + dx / 2;
        const T yint = fint(xhi) - fint(xlo);

Replacement in CarpetX/src/prolongate_3d_rf2.cxx at line 345 [4.1186]

B:BD[3.3518] → [3.3518:4210]

B:BD[3.4210] → [2.133:177]

∅:D[2.177] → [3.4249:5008]

B:BD[3.4249] → [3.4249:5008]

        array<T, n + 3> xs, ys;
        xs[0] = xs[n + 2] = 0 / T(0);
        ys[0] = ys[n + 2] = 0 / T(0);
        const int i0 = n / 2;
        for (int i = -1; i < n; ++i) {
          T x = (i - i0) + CENTERING / T(2);
          T y = f(x);
          xs[i + 2] = x;
          ys[i + 2] = y;
        }
        array<T, 3> x1;
        array<T, 3> y1;
        for (int off = -1; off < 2; ++off) {
          x1[off + 1] = CENTERING / T(4) + off / T(2);
          if (off < 0)
            y1[off + 1] =
                interp1d<CENTERING, CONS, ORDER>()(&ys[i0 + 1], 1, off + 2);
          else
            y1[off + 1] =
                interp1d<CENTERING, CONS, ORDER>()(&ys[i0 + 2], 1, off);
          assert(!CCTK_isnan(y1[off + 1]));
        }
        T dx = x1[1] - x1[0];
        T xlo = x1[0];
        T xhi = x1[2];
        T yint = fint(xhi) - fint(xlo);
        if (!(y1[0] / 4 + y1[1] / 2 + y1[2] / 4 == ys[i0 + 2]) ||
            !(y1[0] * dx / 2 + y1[1] * dx + y1[2] * dx / 2 == yint)) {
          cerr << "settings: CENTERING=" << CENTERING << " ORDER=" << ORDER
               << " order=" << order << "\n";
          cerr << "input:\n";
          for (int i = -1; i < n; ++i)
            cerr << "  xs=" << xs[i + 2] << " ys=" << ys[i + 2] << "\n";
          cerr << "output:\n";
          for (int off = -1; off < 2; ++off)
            cerr << "  x1=" << x1[off + 1] << " y1=" << y1[off + 1] << "\n";
          cerr << "xlo=" << xlo << " xhi=" << xhi << " yint=" << yint << "\n";

[3.3518]

[3.5008]

        // Don't test this, the case (VC,CONS) should not be used
        if (false) {
          array<T, n + 3> xs, ys;
          xs[0] = xs[n + 2] = 0 / T(0);
          ys[0] = ys[n + 2] = 0 / T(0);
          const int i0 = n / 2;
          for (int i = -1; i < n; ++i) {
            const T x = (i - i0) + CENTERING / T(2);
            // T y = f(x);
            const T dx = 1;
            const T xlo = x - dx / 2;
            const T xhi = x + dx / 2;
            const T y = fint(xhi) - fint(xlo);
            xs[i + 2] = x;
            ys[i + 2] = y;
          }
          array<T, 3> x1;
          array<T, 3> y1;
          for (int off = -1; off < 2; ++off) {
            x1[off + 1] = CENTERING / T(4) + off / T(2);
            if (off < 0)
              y1[off + 1] =
                  interp1d<CENTERING, CONS, ORDER>()(&ys[i0 + 1], 1, off + 2);
            else
              y1[off + 1] =
                  interp1d<CENTERING, CONS, ORDER>()(&ys[i0 + 2], 1, off);
            assert(!CCTK_isnan(y1[off + 1]));
          }
          const T dx = x1[1] - x1[0];
          const T xlo = x1[0];
          const T xhi = x1[2];
          const T yint = fint(xhi) - fint(xlo);
          if (!(y1[0] / 4 + y1[1] / 2 + y1[2] / 4 == ys[i0 + 2]) ||
              !(y1[0] * dx / 2 + y1[1] * dx + y1[2] * dx / 2 == yint)) {
            cerr << "settings: CENTERING=" << CENTERING << " ORDER=" << ORDER
                 << " order=" << order << "\n";
            cerr << "input:\n";
            for (int i = -1; i < n; ++i)
              cerr << "  xs=" << xs[i + 2] << " ys=" << ys[i + 2] << "\n";
            cerr << "output:\n";
            for (int off = -1; off < 2; ++off)
              cerr << "  x1=" << x1[off + 1] << " y1=" << y1[off + 1] << "\n";
            cerr << "xlo=" << xlo << " xhi=" << xhi << " yint=" << yint << "\n";
          }
          assert(y1[0] / 4 + y1[1] / 2 + y1[2] / 4 == ys[i0 + 2]);
          assert(y1[0] * dx / 2 + y1[1] * dx + y1[2] * dx / 2 == yint);

Deletion in CarpetX/src/prolongate_3d_rf2.cxx at line 392 [4.1186]

B:BD[3.5018] → [3.5018:5153]

        assert(y1[0] / 4 + y1[1] / 2 + y1[2] / 4 == ys[i0 + 2]);
        assert(y1[0] * dx / 2 + y1[1] * dx + y1[2] * dx / 2 == yint);