Add fdot_abs_integrated_moms diagnostic

core/unit/ctest_array_reduce.c

@@ -102,9 +102,59 @@ void test_sum_reduce_range()
   gkyl_array_release(arr);
 }
 
+void test_sum_abs_reduce_range()
+{
+  int shape[] = {2, 3};
+  struct gkyl_range range;
+  gkyl_range_init_from_shape(&range, 2, shape);
+
+  struct gkyl_array *arr = gkyl_array_new(GKYL_DOUBLE, 2, range.volume);
+  for (size_t i=0; i<arr->size; ++i) {
+    double *d = gkyl_array_fetch(arr, i);
+    d[0] = i%2 == 0 ? 0.5 : -0.5;
+    d[1] = i%2 == 0 ? -1.5 : 1.5;
+  }
+
+  double asum[2];
+  gkyl_array_reduce_range(asum, arr, GKYL_SUM_ABS, &range);
+
+  TEST_CHECK( asum[0] == 0.5*range.volume );
+  TEST_CHECK( asum[1] == 1.5*range.volume );
+
+  gkyl_array_release(arr);
+}
+
 // CUDA specific tests
 #ifdef GKYL_HAVE_CUDA
 
+void test_cu_sum_abs_reduce_range()
+{
+  int shape[] = {2, 3};
+  struct gkyl_range range;
+  gkyl_range_init_from_shape(&range, 2, shape);
+
+  struct gkyl_array *arr = gkyl_array_new(GKYL_DOUBLE, 2, range.volume);
+  struct gkyl_array *arr_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, 2, range.volume);
+  for (size_t i=0; i<arr->size; ++i) {
+    double *d = gkyl_array_fetch(arr, i);
+    d[0] = i%2 == 0 ? 0.5 : -0.5;
+    d[1] = i%2 == 0 ? -1.5 : 1.5;
+  }
+  gkyl_array_copy(arr_cu, arr);
+
+  double asum[2];
+  double *asum_cu = gkyl_cu_malloc(2*sizeof(double));
+  gkyl_array_reduce_range(asum_cu, arr_cu, GKYL_SUM_ABS, &range);
+  gkyl_cu_memcpy(asum, asum_cu, sizeof(asum), GKYL_CU_MEMCPY_D2H);
+
+  TEST_CHECK( asum[0] == 0.5*range.volume );
+  TEST_CHECK( asum[1] == 1.5*range.volume );
+
+  gkyl_cu_free(asum_cu);
+  gkyl_array_release(arr_cu);
+  gkyl_array_release(arr);
+}
+
 void test_cu_array_reduce_max()
 {
 
@@ -374,7 +424,9 @@ TEST_LIST = {
   { "array_reduce", test_reduce },
   { "array_reduce_range", test_reduce_range },
   { "array_reduce_sum_range", test_sum_reduce_range },
+  { "array_reduce_sum_abs_range", test_sum_abs_reduce_range },
 #ifdef GKYL_HAVE_CUDA
+  { "cu_array_reduce_sum_abs_range", test_cu_sum_abs_reduce_range },
   { "cu_array_reduce_max", test_cu_array_reduce_max },
   { "cu_array_reduce_max_big", test_cu_array_reduce_max_big },
   { "cu_array_reduce_range_1d_max", test_cu_array_reduce_range_1d_max },

core/zero/array_dg_reduce.c

@@ -21,6 +21,9 @@ gkyl_array_dg_reducec(double *out, const struct gkyl_array *arr, int comp,
       case GKYL_SUM:                                      
         gkyl_array_dg_reducec_sum_cu(out, arr, comp, basis);
         break;
+      default:
+        assert(false);
+        break;
     }
     return;
   }
@@ -66,6 +69,9 @@ gkyl_array_dg_reducec(double *out, const struct gkyl_array *arr, int comp,
         }
       }
       break;
+    default:
+      assert(false);
+      break;
   }
 }
 
@@ -87,6 +93,9 @@ gkyl_array_dg_reducec_range(double *out, const struct gkyl_array *arr, int comp,
       case GKYL_SUM:
         gkyl_array_dg_reducec_range_sum_cu(out, arr, comp, basis, range);
         break;
+      default:
+        assert(false);
+        break;
     }
     return;
   }
@@ -138,6 +147,8 @@ gkyl_array_dg_reducec_range(double *out, const struct gkyl_array *arr, int comp,
         }
       }
       break;
+    default:
+      assert(false);
+      break;
   }
 }
-

core/zero/array_reduce.c

@@ -20,6 +20,9 @@ gkyl_array_reduce(double *out, const struct gkyl_array *arr, enum gkyl_array_op
       case GKYL_SUM:
         gkyl_array_reduce_sum_cu(out, arr);
         break;
+      default:
+        assert(false);
+        break;
     }
     return;
   }
@@ -55,6 +58,9 @@ gkyl_array_reduce(double *out, const struct gkyl_array *arr, enum gkyl_array_op
           out[k] += d[k];
       }
       break;
+    default:
+      assert(false);
+      break;
   }
 }
 
@@ -76,6 +82,9 @@ gkyl_array_reduce_range(double *res,
       case GKYL_SUM:
         gkyl_array_reduce_range_sum_cu(res, arr, range);
         break;
+      case GKYL_SUM_ABS:
+        gkyl_array_reduce_range_sum_abs_cu(res, arr, range);
+        break;
     }
     return;
   }
@@ -116,6 +125,15 @@ gkyl_array_reduce_range(double *res,
           res[i] += d[i];
       }
       break;
+    case GKYL_SUM_ABS:
+      for (long i=0; i<n; ++i) res[i] = 0;
+
+      while (gkyl_range_iter_next(&iter)) {
+        long start = gkyl_range_idx(range, iter.idx);
+        const double *d = gkyl_array_cfetch(arr, start);
+        for (long i=0; i<n; ++i)
+          res[i] += fabs(d[i]);
+      }
+      break;
   }
 }
-

core/zero/array_reduce_cu.cu

@@ -327,3 +327,49 @@ gkyl_array_reduce_range_sum_cu(double *out_d, const struct gkyl_array* inp, cons
   cudaDeviceSynchronize();
 }
 
+template <unsigned int BLOCKSIZE>
+__global__ void
+arraySumAbs_range_blockRedAtomic_cub(const struct gkyl_array* inp,
+  const struct gkyl_range range, double* out)
+{
+  unsigned long linc = blockIdx.x*blockDim.x + threadIdx.x;
+
+  typedef cub::BlockReduce<double, BLOCKSIZE> BlockReduceT;
+  __shared__ typename BlockReduceT::TempStorage temp;
+
+  long nCells = range.volume;
+  size_t nComp = inp->ncomp;
+  int idx[GKYL_MAX_DIM];
+
+  for (size_t k = 0; k < nComp; ++k) {
+    double f = 0;
+    if (linc < nCells) {
+      gkyl_sub_range_inv_idx(&range, linc, idx);
+      long start = gkyl_range_idx(&range, idx);
+      const double* fptr = (const double*) gkyl_array_cfetch(inp, start);
+      f = fabs(fptr[k]);
+    }
+    double bResult = BlockReduceT(temp).Reduce(f,
+#if CUDART_VERSION > 12090
+      ::cuda::std::plus()
+#else
+      cub::Sum()
+#endif
+    );
+    if (threadIdx.x == 0)
+      atomicAdd(&out[k], bResult);
+  }
+}
+
+void
+gkyl_array_reduce_range_sum_abs_cu(double *out_d, const struct gkyl_array* inp,
+  const struct gkyl_range *range)
+{
+  gkyl_cu_memset(out_d, 0, inp->ncomp*sizeof(double));
+
+  const int nthreads = GKYL_DEFAULT_NUM_THREADS;
+  int nblocks = gkyl_int_div_up(range->volume, nthreads);
+  arraySumAbs_range_blockRedAtomic_cub<nthreads><<<nblocks, nthreads>>>(inp->on_dev, *range, out_d);
+  // Device synchronize required because out_d may be host pinned memory.
+  cudaDeviceSynchronize();
+}

core/zero/gkyl_array_reduce.h

@@ -24,4 +24,3 @@ void gkyl_array_reduce(double *res, const struct gkyl_array *arr, enum gkyl_arra
  */
 void gkyl_array_reduce_range(double *res,
   const struct gkyl_array *arr, enum gkyl_array_op op, const struct gkyl_range *range);
-

core/zero/gkyl_array_reduce_priv.h

@@ -55,4 +55,14 @@ void gkyl_array_reduce_sum_cu(double *out_d, const struct gkyl_array* inp);
  */
 void gkyl_array_reduce_range_sum_cu(double *out_d, const struct gkyl_array* inp, const struct gkyl_range *range);
 
+/**
+ * Sum the absolute values of a gkyl_array component-wise over a specified range.
+ *
+ * @param out_d A device array with as many elements as the 'inp' has components.
+ * @param inp A gkyl_array to be reduced.
+ * @param range A gkyl_range over which to perform the reduction.
+ */
+void gkyl_array_reduce_range_sum_abs_cu(double *out_d, const struct gkyl_array* inp,
+  const struct gkyl_range *range);
+
 #endif

core/zero/gkyl_elem_type.h

@@ -7,7 +7,7 @@
 enum gkyl_elem_type { GKYL_INT, GKYL_FLOAT, GKYL_DOUBLE, GKYL_INT_64, GKYL_LONG, GKYL_USER };
 
 // Array reduce operators
-enum gkyl_array_op { GKYL_MIN, GKYL_MAX, GKYL_SUM };
+enum gkyl_array_op { GKYL_MIN, GKYL_MAX, GKYL_SUM, GKYL_SUM_ABS };
 
 // file types for raw IO of Gkeyll data
 enum gkyl_file_type {

gkeyll/lua/DataStruct/ZeroArray.lua

@@ -95,7 +95,7 @@ struct gkyl_array* gkyl_array_cu_host_new(enum gkyl_elem_type type, size_t ncomp
  */
 bool gkyl_array_is_cu_dev(const struct gkyl_array *arr);
 
-enum gkyl_array_op { GKYL_MIN, GKYL_MAX, GKYL_SUM };
+enum gkyl_array_op { GKYL_MIN, GKYL_MAX, GKYL_SUM, GKYL_SUM_ABS };
 
 
 /**
@@ -519,6 +519,8 @@ local array_fn = {
          enum = 1
       elseif op == "sum" then
          enum = 2
+      elseif op == "sum_abs" then
+         enum = 3
       end
       ffiC.gkyl_array_reduce_range(out, self, enum, rng)
    end,