Add/clean core and gyrokinetic unit tests

core/unit/ctest_array_ops.c

@@ -1142,6 +1142,118 @@ void test_array_min_range_ho() {
   test_array_min_range(false);
 }
 
+void
+test_array_new_meta()
+{
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 3, 10);
+  TEST_CHECK( a->type == GKYL_DOUBLE );
+  TEST_CHECK( a->ncomp == 3 );
+  TEST_CHECK( a->size == 10 );
+  TEST_CHECK( gkyl_array_is_cu_dev(a) == false );
+  gkyl_array_release(a);
+}
+
+void
+test_array_clear_basic()
+{
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 1, 50);
+  gkyl_array_clear(a, 3.5);
+  double *d = a->data;
+  for (unsigned i=0; i<a->size*a->ncomp; ++i)
+    TEST_CHECK( d[i] == 3.5 );
+  gkyl_array_release(a);
+}
+
+void
+test_array_scale_basic()
+{
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 1, 20);
+  gkyl_array_clear(a, 2.0);
+  gkyl_array_scale(a, 4.0);
+  double *d = a->data;
+  for (unsigned i=0; i<a->size; ++i)
+    TEST_CHECK( gkyl_compare_double(d[i], 8.0, 1e-14) );
+  gkyl_array_release(a);
+}
+
+void
+test_array_accumulate_basic()
+{
+  // out = out + a*inp
+  struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 1, 20);
+  struct gkyl_array *inp = gkyl_array_new(GKYL_DOUBLE, 1, 20);
+  gkyl_array_clear(out, 1.0);
+  gkyl_array_clear(inp, 2.0);
+  gkyl_array_accumulate(out, 3.0, inp);
+  double *d = out->data;
+  for (unsigned i=0; i<out->size; ++i)
+    TEST_CHECK( gkyl_compare_double(d[i], 1.0 + 3.0*2.0, 1e-14) );
+  gkyl_array_release(out); gkyl_array_release(inp);
+}
+
+void
+test_array_set_basic()
+{
+  // out = a*inp
+  struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 1, 20);
+  struct gkyl_array *inp = gkyl_array_new(GKYL_DOUBLE, 1, 20);
+  gkyl_array_clear(out, 99.0);
+  gkyl_array_clear(inp, 5.0);
+  gkyl_array_set(out, 2.0, inp);
+  double *d = out->data;
+  for (unsigned i=0; i<out->size; ++i)
+    TEST_CHECK( gkyl_compare_double(d[i], 10.0, 1e-14) );
+  gkyl_array_release(out); gkyl_array_release(inp);
+}
+
+void
+test_array_shiftc_basic()
+{
+  // shift component k by a (out[k] += a)
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 3, 10);
+  gkyl_array_clear(a, 0.0);
+  gkyl_array_shiftc(a, 7.0, 1); // shift component 1
+  for (unsigned i=0; i<a->size; ++i) {
+    double *row = gkyl_array_fetch(a, i);
+    TEST_CHECK( row[0] == 0.0 );
+    TEST_CHECK( gkyl_compare_double(row[1], 7.0, 1e-14) );
+    TEST_CHECK( row[2] == 0.0 );
+  }
+  gkyl_array_release(a);
+}
+
+void
+test_array_copy_basic()
+{
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 2, 15);
+  struct gkyl_array *b = gkyl_array_new(GKYL_DOUBLE, 2, 15);
+  double *ad = a->data;
+  for (unsigned i=0; i<a->size*a->ncomp; ++i) ad[i] = 0.5*i;
+  gkyl_array_copy(b, a);
+  double *bd = b->data;
+  for (unsigned i=0; i<b->size*b->ncomp; ++i)
+    TEST_CHECK( bd[i] == 0.5*i );
+  gkyl_array_release(a); gkyl_array_release(b);
+}
+
+void
+test_array_fetch_multicomp()
+{
+  // Write per-component values via fetch and read back.
+  struct gkyl_array *a = gkyl_array_new(GKYL_DOUBLE, 4, 8);
+  for (unsigned i=0; i<a->size; ++i) {
+    double *row = gkyl_array_fetch(a, i);
+    for (unsigned c=0; c<a->ncomp; ++c)
+      row[c] = 100.0*i + c;
+  }
+  for (unsigned i=0; i<a->size; ++i) {
+    const double *row = gkyl_array_cfetch(a, i);
+    for (unsigned c=0; c<a->ncomp; ++c)
+      TEST_CHECK( row[c] == 100.0*i + c );
+  }
+  gkyl_array_release(a);
+}
+
 // Cuda specific tests
 #ifdef GKYL_HAVE_CUDA
 
@@ -2184,6 +2296,14 @@ TEST_LIST = {
   { "array_copy_range", test_array_copy_range},
   { "array_copy_split", test_array_copy_split },
   { "array_copy_range_to_range_diff_range_dim", test_array_copy_range_to_range_diff_range_dim_ho},
+  { "array_new_meta", test_array_new_meta },
+  { "array_clear_basic", test_array_clear_basic },
+  { "array_scale_basic", test_array_scale_basic },
+  { "array_accumulate_basic", test_array_accumulate_basic },
+  { "array_set_basic", test_array_set_basic },
+  { "array_shiftc_basic", test_array_shiftc_basic },
+  { "array_copy_basic", test_array_copy_basic },
+  { "array_fetch_multicomp", test_array_fetch_multicomp },
 #ifdef GKYL_HAVE_CUDA
   { "cu_array_clear", test_cu_array_clear},
   { "cu_array_clear_range", test_cu_array_clear_range},

core/unit/ctest_basis.c

@@ -8,6 +8,130 @@
 #include <gkyl_basis_gkhyb_1x2v_p1_surfx3_eval_quad.h>
 #include <gkyl_basis_gkhyb_1x2v_p1_upwind_quad_to_modal.h>
 
+static int ipow(int b, int e) { int r=1; for (int i=0;i<e;++i) r*=b; return r; }
+
+void
+test_basis_num_basis_serendip()
+{
+  struct gkyl_basis b;
+
+  // 1D serendipity: num_basis = poly_order+1
+  gkyl_cart_modal_serendip(&b, 1, 1); TEST_CHECK( b.num_basis == 2 );
+  gkyl_cart_modal_serendip(&b, 1, 2); TEST_CHECK( b.num_basis == 3 );
+  gkyl_cart_modal_serendip(&b, 1, 3); TEST_CHECK( b.num_basis == 4 );
+
+  // 2D serendipity: p1=4, p2=8, p3=12
+  gkyl_cart_modal_serendip(&b, 2, 1); TEST_CHECK( b.num_basis == 4 );
+  gkyl_cart_modal_serendip(&b, 2, 2); TEST_CHECK( b.num_basis == 8 );
+  gkyl_cart_modal_serendip(&b, 2, 3); TEST_CHECK( b.num_basis == 12 );
+
+  // 3D serendipity: p1=8, p2=20, p3=32
+  gkyl_cart_modal_serendip(&b, 3, 1); TEST_CHECK( b.num_basis == 8 );
+  gkyl_cart_modal_serendip(&b, 3, 2); TEST_CHECK( b.num_basis == 20 );
+
+  // metadata sanity
+  gkyl_cart_modal_serendip(&b, 2, 2);
+  TEST_CHECK( b.ndim == 2 );
+  TEST_CHECK( b.poly_order == 2 );
+  TEST_CHECK( b.b_type == GKYL_BASIS_MODAL_SERENDIPITY );
+}
+
+void
+test_basis_num_basis_tensor()
+{
+  struct gkyl_basis b;
+  // tensor basis: num_basis = (poly_order+1)^ndim
+  for (int ndim=1; ndim<=3; ++ndim) {
+    for (int p=1; p<=2; ++p) {
+      gkyl_cart_modal_tensor(&b, ndim, p);
+      TEST_CHECK( b.num_basis == (unsigned) ipow(p+1, ndim) );
+      TEST_CHECK( b.b_type == GKYL_BASIS_MODAL_TENSOR );
+    }
+  }
+}
+
+void
+test_basis_eval_expand_consistency()
+{
+  // eval_expand(z,f) must equal sum_i f[i]*b_i(z), where b = eval(z).
+  struct gkyl_basis b;
+  gkyl_cart_modal_serendip(&b, 2, 2);
+
+  double z[2] = { 0.3, -0.7 };
+  double bz[8];
+  b.eval(z, bz);
+
+  double f[8] = { 1.0, 2.0, -1.5, 0.5, 3.0, -2.0, 0.25, -0.75 };
+  double ev = b.eval_expand(z, f);
+
+  double sum = 0.0;
+  for (int i=0; i<8; ++i) sum += f[i]*bz[i];
+  TEST_CHECK( gkyl_compare_double(ev, sum, 1e-13) );
+}
+
+void
+test_basis_const_mode()
+{
+  // The first (0th) basis function is constant: same value at any z, and >0.
+  struct gkyl_basis b;
+  gkyl_cart_modal_serendip(&b, 1, 2);
+
+  double z1[1] = { -0.9 }, z2[1] = { 0.4 };
+  double b1[3], b2[3];
+  b.eval(z1, b1);
+  b.eval(z2, b2);
+  TEST_CHECK( gkyl_compare_double(b1[0], b2[0], 1e-14) );
+  TEST_CHECK( b1[0] > 0.0 );
+
+  // Pure constant expansion evaluates to the same value everywhere.
+  double f[3] = { 1.0, 0.0, 0.0 };
+  TEST_CHECK( gkyl_compare_double(b.eval_expand(z1, f), b.eval_expand(z2, f), 1e-14) );
+}
+
+void
+test_basis_flip_odd_involution()
+{
+  // Applying flip_odd_sign twice in the same direction recovers the input.
+  struct gkyl_basis b;
+  gkyl_cart_modal_serendip(&b, 2, 2);
+
+  double f[8] = { 1.0, 2.0, -1.5, 0.5, 3.0, -2.0, 0.25, -0.75 };
+  double f1[8], f2[8];
+  b.flip_odd_sign(0, f, f1);
+  b.flip_odd_sign(0, f1, f2);
+  for (int i=0; i<8; ++i)
+    TEST_CHECK( gkyl_compare_double(f2[i], f[i], 1e-14) );
+}
+
+void
+test_basis_flip_odd_reflection()
+{
+  // flip_odd_sign(dir) on coefficients corresponds to reflecting z->-z in dir:
+  // expand(flip(f))(z) == expand(f)(z with z[dir] negated).
+  struct gkyl_basis b;
+  gkyl_cart_modal_serendip(&b, 1, 2);
+
+  double f[3] = { 1.5, -0.5, 0.8 };
+  double ff[3];
+  b.flip_odd_sign(0, f, ff);
+
+  double z[1] = { 0.6 }, zm[1] = { -0.6 };
+  TEST_CHECK( gkyl_compare_double(b.eval_expand(z, ff), b.eval_expand(zm, f), 1e-13) );
+}
+
+void
+test_basis_new_release()
+{
+  struct gkyl_basis *b = gkyl_cart_modal_serendip_new(3, 1);
+  TEST_CHECK( b->num_basis == 8 );
+  TEST_CHECK( b->ndim == 3 );
+  gkyl_cart_modal_basis_release(b);
+
+  struct gkyl_basis *t = gkyl_cart_modal_tensor_new(2, 2);
+  TEST_CHECK( t->num_basis == 9 );
+  gkyl_cart_modal_basis_release(t);
+}
+
 void
 test_ser_1d_p0_members(struct gkyl_basis basis1)
 {
@@ -753,6 +877,13 @@ TEST_LIST = {
   { "ten_2d", test_ten_2d },
   { "hyb", test_hyb },
   { "gkhyb", test_gkhyb },
+  { "num_basis_serendip", test_basis_num_basis_serendip },
+  { "num_basis_tensor", test_basis_num_basis_tensor },
+  { "eval_expand_consistency", test_basis_eval_expand_consistency },
+  { "const_mode", test_basis_const_mode },
+  { "flip_odd_involution", test_basis_flip_odd_involution },
+  { "flip_odd_reflection", test_basis_flip_odd_reflection },
+  { "basis_new_release", test_basis_new_release },
 #ifdef GKYL_HAVE_CUDA
   { "cu_ser_2d", test_cu_ser_2d },
 #endif    

core/unit/ctest_eval_on_nodes.c

@@ -275,7 +275,7 @@ void test_1x1v_hyb(int poly_order, int test_func_op)
 
   // project distribution function on basis
   gkyl_eval_on_nodes_advance(evup, 0.0, &arr_range, distf);
-  gkyl_grid_sub_array_write(&grid, &arr_range, 0, distf, "ctest_eval_on_nodes_distf_hyb.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &arr_range, 0, distf, "ctest_eval_on_nodes_distf_hyb.gkyl");
 
   double *dfll = gkyl_array_fetch(distf, 0);  // left, low cell
   double *dflu = gkyl_array_fetch(distf, 1);  // left, up cell
@@ -1149,7 +1149,7 @@ void test_2x2v_hyb(int poly_order, int test_func_op)
 
   // project distribution function on basis
   gkyl_eval_on_nodes_advance(evup, 0.0, &arr_range, distf);
-  gkyl_grid_sub_array_write(&grid, &arr_range, 0, distf, "ctest_eval_on_nodes_distf_hyb.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &arr_range, 0, distf, "ctest_eval_on_nodes_distf_hyb.gkyl");
 
   if (poly_order == 1) {
     if (test_func_op==0) {

core/unit/ctest_gauss_quad_data.c

@@ -0,0 +1,99 @@
+// Unit tests for Gauss-Legendre quadrature data in gkyl_gauss_quad_data.h
+#include <acutest.h>
+#include <gkyl_gauss_quad_data.h>
+#include <gkyl_util.h>
+#include <math.h>
+
+// Integrate x^p over [-1,1] using N-point Gauss-Legendre quadrature.
+static double
+quad_int_monomial(int N, int p)
+{
+  const double *x = gkyl_gauss_ordinates[N];
+  const double *w = gkyl_gauss_weights[N];
+  double sum = 0.0;
+  for (int i=0; i<N; ++i)
+    sum += w[i]*pow(x[i], p);
+  return sum;
+}
+
+void
+test_gauss_weights_sum()
+{
+  // Sum of weights = integral of 1 over [-1,1] = 2, for every N.
+  for (int N=1; N<=gkyl_gauss_max; ++N) {
+    const double *w = gkyl_gauss_weights[N];
+    double sum = 0.0;
+    for (int i=0; i<N; ++i) sum += w[i];
+    TEST_CHECK( gkyl_compare_double(sum, 2.0, 1e-13) );
+    TEST_MSG("N=%d sum=%g", N, sum);
+  }
+}
+
+void
+test_gauss_ordinates_symmetric()
+{
+  // Ordinates are symmetric about 0; reversed list = negated list.
+  for (int N=1; N<=gkyl_gauss_max; ++N) {
+    const double *x = gkyl_gauss_ordinates[N];
+    for (int i=0; i<N; ++i)
+      TEST_CHECK( gkyl_compare_double(x[i], -x[N-1-i], 1e-13) );
+  }
+}
+
+void
+test_gauss_ordinates_in_range()
+{
+  // All ordinates lie strictly inside (-1,1).
+  for (int N=1; N<=gkyl_gauss_max; ++N) {
+    const double *x = gkyl_gauss_ordinates[N];
+    for (int i=0; i<N; ++i)
+      TEST_CHECK( x[i] > -1.0 && x[i] < 1.0 );
+  }
+}
+
+void
+test_gauss_exactness()
+{
+  // N-point Gauss-Legendre integrates polynomials up to degree 2N-1 exactly.
+  // Integral of x^p over [-1,1] = 0 (odd p) or 2/(p+1) (even p).
+  for (int N=1; N<=gkyl_gauss_max; ++N) {
+    int maxdeg = 2*N - 1;
+    for (int p=0; p<=maxdeg; ++p) {
+      double exact = (p%2==1) ? 0.0 : 2.0/(p+1);
+      double approx = quad_int_monomial(N, p);
+      TEST_CHECK( gkyl_compare_double(approx, exact, 1e-11) );
+      TEST_MSG("N=%d p=%d exact=%g approx=%g", N, p, exact, approx);
+    }
+  }
+}
+
+void
+test_gauss_specific_values()
+{
+  // 2-point Gauss: ordinates +/- 1/sqrt(3), weights 1,1.
+  const double *x2 = gkyl_gauss_ordinates[2];
+  const double *w2 = gkyl_gauss_weights[2];
+  TEST_CHECK( gkyl_compare_double(fabs(x2[0]), 1.0/sqrt(3.0), 1e-12) );
+  TEST_CHECK( gkyl_compare_double(w2[0], 1.0, 1e-12) );
+  TEST_CHECK( gkyl_compare_double(w2[1], 1.0, 1e-12) );
+
+  // 1-point Gauss: ordinate 0, weight 2.
+  TEST_CHECK( gkyl_compare_double(gkyl_gauss_ordinates[1][0], 0.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_gauss_weights[1][0], 2.0, 1e-13) );
+
+  // 3-point Gauss: middle ordinate 0, weights 5/9, 8/9, 5/9.
+  const double *x3 = gkyl_gauss_ordinates[3];
+  const double *w3 = gkyl_gauss_weights[3];
+  TEST_CHECK( gkyl_compare_double(x3[1], 0.0, 1e-12) );
+  TEST_CHECK( gkyl_compare_double(w3[1], 8.0/9.0, 1e-12) );
+  TEST_CHECK( gkyl_compare_double(w3[0], 5.0/9.0, 1e-12) );
+}
+
+TEST_LIST = {
+  { "gauss_weights_sum", test_gauss_weights_sum },
+  { "gauss_ordinates_symmetric", test_gauss_ordinates_symmetric },
+  { "gauss_ordinates_in_range", test_gauss_ordinates_in_range },
+  { "gauss_exactness", test_gauss_exactness },
+  { "gauss_specific_values", test_gauss_specific_values },
+  { NULL, NULL },
+};