C++ 8th Order Divergence

src/cpp/divergence.cpp

@@ -18,7 +18,7 @@
 // 1-D Constructor
 Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
   assert(!(k % 2));
-  assert(k > 1 && k < 7);
+  assert(k > 1 && k < 9);
   assert(m > 2 * k);
 
   switch (k) {
@@ -100,6 +100,88 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       , 1706.0 / 1457.0 , 3124.0 / 5901.0 , 887.0 / 531.0 , 929.0 / 2002.0
       , 2383.0 / 2005.0 };
     break;
+  case 8:
+    /*
+    WARNING:
+      At the 8th order, the weight matrix Q loses positive definiteness, 
+      so the inner product induced by Q is no longer well-defined. If 
+      the inner product is not valid, the discrete integration by parts 
+      identity has no meaning, which breaks the structure that makes 
+      the divergence mimetic.
+    */
+    // A
+    at(1, 0) = -1423.0 / 1792.0;
+    at(1, 1) = -491.0 / 7168.0;
+    at(1, 2) = 7753.0 / 3072.0;
+    at(1, 3) = -18509.0 / 5120.0;
+    at(1, 4) = 3535.0 / 1024.0;
+    at(1, 5) = -2279.0 / 1024.0;
+    at(1, 6) = 953.0 / 1024.0;
+    at(1, 7) = -1637.0 / 7168.0;
+    at(1, 8) = 2689.0 / 107520.0;
+    at(2, 0) = 2689.0 / 107520.0;
+    at(2, 1) = -36527.0 / 35840.0;
+    at(2, 2) = 4259.0 / 5120.0;
+    at(2, 3) = 6497.0 / 15360.0;
+    at(2, 4) = -475.0 / 1024.0;
+    at(2, 5) = 1541.0 / 5120.0;
+    at(2, 6) = -639.0 / 5120.0;
+    at(2, 7) = 1087.0 / 35840.0;
+    at(2, 8) = -59.0 / 17920.0;
+    at(3, 0) = -59.0 / 17920.0;
+    at(3, 1) = 1175.0 / 21504.0;
+    at(3, 2) = -1165.0 / 1024.0;
+    at(3, 3) = 1135.0 / 1024.0;
+    at(3, 4) = 25.0 / 3072.0;
+    at(3, 5) = -251.0 / 5120.0;
+    at(3, 6) = 25.0 / 1024.0;;
+    at(3, 7) = -45.0 / 7168.0;
+    at(3, 8) = 5.0 / 7168.0;
+    // A'
+    at(m, m) = 1423.0 / 1792.0;
+    at(m, m - 1) = 491.0 / 7168.0;
+    at(m, m - 2) = -7753.0 / 3072.0;
+    at(m, m - 3) = 18509.0 / 5120.0;
+    at(m, m - 4) = -3535.0 / 1024.0;
+    at(m, m - 5) = 2279.0 / 1024.0;
+    at(m, m - 6) = -953.0 / 1024.0;
+    at(m, m - 7) = 1637.0 / 7168.0;
+    at(m, m - 8) = -2689.0 / 107520.0;
+    at(m - 1, m) = -2689.0 / 107520.0;
+    at(m - 1, m - 1) = 36527.0 / 35840.0;
+    at(m - 1, m - 2) = -4259.0 / 5120.0;
+    at(m - 1, m - 3) = -6497.0 / 15360.0;
+    at(m - 1, m - 4) = 475.0 / 1024.0;
+    at(m - 1, m - 5) = -1541.0 / 5120.0;
+    at(m - 1, m - 6) = 639.0 / 5120.0;
+    at(m - 1, m - 7) = -1087.0 / 35840.0;
+    at(m - 1, m - 8) = 59.0 / 17920.0;
+    at(m - 2, m) = 59.0 / 17920.0;
+    at(m - 2, m - 1) = -1175.0 / 21504.0;
+    at(m - 2, m - 2) = 1165.0 / 1024.0;
+    at(m - 2, m - 3) = -1135.0 / 1024.0;
+    at(m - 2, m - 4) = -25.0 / 3072.0;
+    at(m - 2, m - 5) = 251.0 / 5120.0;
+    at(m - 2, m - 6) = -25.0 / 1024.0;;
+    at(m - 2, m - 7) = 45.0 / 7168.0;
+    at(m - 2, m - 8) = -5.0 / 7168.0;
+    // Middle
+    for (u32 i = 4; i < m - 2; ++i) {
+      at(i, i - 4) = 5.0 / 7168.0;
+      at(i, i - 3) = -49.0 / 5120.0;
+      at(i, i - 2) = 245.0 / 3072.0;
+      at(i, i - 1) = -1225.0 / 1024.0;
+      at(i, i) = 1225.0 / 1024.0;
+      at(i, i + 1) = -245.0 / 3072.0;
+      at(i, i + 2) = 49.0 / 5120.0;
+      at(i, i + 3) = -5.0 / 7168.0;
+    }
+    // Weights
+    Q  = { 1558.0 / 1247.0 , 271.0 / 3660.0 , 3225.0 / 1181.0 , -1103.0 / 1050.0
+      , 797.0 / 312.0 , 632.0 / 2273.0 , 755.0 / 641.0 , 859.0 / 869.0
+      , 966.0 / 971.0 , 859.0 / 869.0 , 755.0 / 641.0 , 632.0 / 2273.0
+      , 797.0 / 312.0 , -1103.0 / 1050.0 , 3225.0 / 1181.0 , 271.0 / 3660.0
+      , 1558.0 / 1247.0 };
   }
 
   // Scaling

src/cpp/divergence.h

@@ -22,6 +22,12 @@
 /**
  * @brief Mimetic Divergence operator
  *
+ * WARNING:
+ *    At the 8th order, the weight matrix Q loses positive definiteness, 
+ *    so the inner product induced by Q is no longer well-defined. If 
+ *    the inner product is not valid, the discrete integration by parts 
+ *    identity has no meaning, which breaks the structure that makes 
+ *    the divergence mimetic.
  */
 class Divergence : public sp_mat {
 

tests/cpp/test1.cpp

@@ -15,7 +15,7 @@ void run_nullity_test(int k, Real tol) {
 
 TEST(DivergenceTests, Nullity) {
     Real tol = 1e-10;
-    for (int k : {2, 4, 6}) {
+    for (int k : {2, 4, 6, 8}) {
         run_nullity_test(k, tol);
     }
 }

tests/cpp/test3.cpp

@@ -15,7 +15,7 @@ void run_nullity_test(int k, Real tol) {
 
 TEST(LaplacianTests, Nullity) {
     Real tol = 1e-10;
-    for (int k : {2, 4, 6}) {
+    for (int k : {2, 4, 6, 8}) {
         run_nullity_test(k, tol);
     }
 }

tests/cpp/test5.cpp

@@ -46,7 +46,7 @@ void run_test(int k, vec grid_sizes) {
 
 TEST(PoissonTests, Accuracy) {
     vec grid_sizes = {20, 40};
-    for (int k : {2, 4, 6}) {
+    for (int k : {2, 4, 6}) { // We ignore k=8 b/c e = O(1e-12)
         run_test(k, grid_sizes);
     }
 }