diff --git a/core/unit/ctest_array_ops.c b/core/unit/ctest_array_ops.c
index 2da77cac69..7a8346a9c4 100644
--- a/core/unit/ctest_array_ops.c
+++ b/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},
diff --git a/core/unit/ctest_basis.c b/core/unit/ctest_basis.c
index 84bf46c787..7faa4ec503 100644
--- a/core/unit/ctest_basis.c
+++ b/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    
diff --git a/core/unit/ctest_eval_on_nodes.c b/core/unit/ctest_eval_on_nodes.c
index 92caea6935..4b275ec53d 100644
--- a/core/unit/ctest_eval_on_nodes.c
+++ b/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) {
diff --git a/core/unit/ctest_gauss_quad_data.c b/core/unit/ctest_gauss_quad_data.c
new file mode 100644
index 0000000000..8f121d62a5
--- /dev/null
+++ b/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 },
+};
diff --git a/core/unit/ctest_mat.c b/core/unit/ctest_mat.c
index 37c42449d7..b14d9b53f9 100644
--- a/core/unit/ctest_mat.c
+++ b/core/unit/ctest_mat.c
@@ -2,6 +2,7 @@
 #include <gkyl_alloc.h>
 #include <gkyl_mat.h>
 #include <gkyl_mat_priv.h>
+#include <gkyl_util.h>
 
 void
 test_mat_base()
@@ -885,14 +886,158 @@ void test_cu_mat_mm_arrays()
 
 #endif
 
+void
+test_mat_new_set_get()
+{
+  struct gkyl_mat *m = gkyl_mat_new(3, 2, 7.0);
+  TEST_CHECK( m->nr == 3 );
+  TEST_CHECK( m->nc == 2 );
+  // initial value
+  for (size_t r=0; r<3; ++r)
+    for (size_t c=0; c<2; ++c)
+      TEST_CHECK( gkyl_mat_get(m, r, c) == 7.0 );
+
+  gkyl_mat_set(m, 0, 0, 1.0);
+  gkyl_mat_set(m, 2, 1, -3.5);
+  TEST_CHECK( gkyl_mat_get(m, 0, 0) == 1.0 );
+  TEST_CHECK( gkyl_mat_get(m, 2, 1) == -3.5 );
+  // unchanged entry
+  TEST_CHECK( gkyl_mat_get(m, 1, 0) == 7.0 );
+
+  gkyl_mat_release(m);
+}
+
+void
+test_mat_clear()
+{
+  struct gkyl_mat *m = gkyl_mat_new(4, 4, 1.0);
+  gkyl_mat_clear(m, 0.0);
+  for (size_t r=0; r<4; ++r)
+    for (size_t c=0; c<4; ++c)
+      TEST_CHECK( gkyl_mat_get(m, r, c) == 0.0 );
+  gkyl_mat_release(m);
+}
+
+void
+test_mat_diag()
+{
+  struct gkyl_mat *m = gkyl_mat_new(3, 3, 9.0);
+  gkyl_mat_diag(m, 2.0);
+  for (size_t r=0; r<3; ++r)
+    for (size_t c=0; c<3; ++c)
+      TEST_CHECK( gkyl_mat_get(m, r, c) == (r==c ? 2.0 : 0.0) );
+  gkyl_mat_release(m);
+}
+
+void
+test_mat_mm_identity()
+{
+  // A * I == A
+  struct gkyl_mat *A = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_set(A, 0, 0, 1.0); gkyl_mat_set(A, 0, 1, 2.0);
+  gkyl_mat_set(A, 1, 0, 3.0); gkyl_mat_set(A, 1, 1, 4.0);
+
+  struct gkyl_mat *I = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_diag(I, 1.0);
+
+  struct gkyl_mat *C = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_mm(1.0, 0.0, GKYL_NO_TRANS, A, GKYL_NO_TRANS, I, C, false);
+
+  for (size_t r=0; r<2; ++r)
+    for (size_t c=0; c<2; ++c)
+      TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, r, c), gkyl_mat_get(A, r, c), 1e-14) );
+
+  gkyl_mat_release(A); gkyl_mat_release(I); gkyl_mat_release(C);
+}
 
+void
+test_mat_mm_known()
+{
+  // [1 2; 3 4] * [5 6; 7 8] = [19 22; 43 50]
+  struct gkyl_mat *A = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_set(A, 0, 0, 1.0); gkyl_mat_set(A, 0, 1, 2.0);
+  gkyl_mat_set(A, 1, 0, 3.0); gkyl_mat_set(A, 1, 1, 4.0);
+
+  struct gkyl_mat *B = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_set(B, 0, 0, 5.0); gkyl_mat_set(B, 0, 1, 6.0);
+  gkyl_mat_set(B, 1, 0, 7.0); gkyl_mat_set(B, 1, 1, 8.0);
+
+  struct gkyl_mat *C = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_mm(1.0, 0.0, GKYL_NO_TRANS, A, GKYL_NO_TRANS, B, C, false);
+
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 0, 0), 19.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 0, 1), 22.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 1, 0), 43.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 1, 1), 50.0, 1e-13) );
+
+  // alpha scaling: 2*(A*B)
+  struct gkyl_mat *C2 = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_mm(2.0, 0.0, GKYL_NO_TRANS, A, GKYL_NO_TRANS, B, C2, false);
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C2, 0, 0), 38.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C2, 1, 1), 100.0, 1e-13) );
+
+  gkyl_mat_release(A); gkyl_mat_release(B);
+  gkyl_mat_release(C); gkyl_mat_release(C2);
+}
+
+void
+test_mat_mm_transpose()
+{
+  // A^T with A = [1 2; 3 4] is [1 3; 2 4]; (A^T)*I checks transpose handling
+  struct gkyl_mat *A = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_set(A, 0, 0, 1.0); gkyl_mat_set(A, 0, 1, 2.0);
+  gkyl_mat_set(A, 1, 0, 3.0); gkyl_mat_set(A, 1, 1, 4.0);
 
+  struct gkyl_mat *I = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_diag(I, 1.0);
+
+  struct gkyl_mat *C = gkyl_mat_new(2, 2, 0.0);
+  gkyl_mat_mm(1.0, 0.0, GKYL_TRANS, A, GKYL_NO_TRANS, I, C, false);
+
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 0, 0), 1.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 0, 1), 3.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 1, 0), 2.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(gkyl_mat_get(C, 1, 1), 4.0, 1e-14) );
+
+  gkyl_mat_release(A); gkyl_mat_release(I); gkyl_mat_release(C);
+}
+
+void
+test_nmat_new()
+{
+  // Batched matrices: 3 matrices each 2x2.
+  struct gkyl_nmat *nm = gkyl_nmat_new(3, 2, 2);
+  TEST_CHECK( nm->num == 3 );
+  TEST_CHECK( nm->nr == 2 );
+  TEST_CHECK( nm->nc == 2 );
+  TEST_CHECK( !gkyl_nmat_is_cu_dev(nm) );
+
+  // Fill each matrix differently and read back.
+  for (size_t k=0; k<3; ++k) {
+    struct gkyl_mat mk = gkyl_nmat_get(nm, k);
+    gkyl_mat_clear(&mk, 0.0);
+    gkyl_mat_set(&mk, 0, 0, (double) k+1);
+  }
+  for (size_t k=0; k<3; ++k) {
+    struct gkyl_mat mk = gkyl_nmat_get(nm, k);
+    TEST_CHECK( gkyl_mat_get(&mk, 0, 0) == (double) k+1 );
+  }
+
+  gkyl_nmat_release(nm);
+}
 
 TEST_LIST = {
   { "mat_base", test_mat_base },
+  { "mat_new_set_get", test_mat_new_set_get },
+  { "mat_clear", test_mat_clear },
+  { "mat_diag", test_mat_diag },
+  { "mat_mm_identity", test_mat_mm_identity },
+  { "mat_mm_known", test_mat_mm_known },
+  { "mat_mm_transpose", test_mat_mm_transpose },
   { "mat_mm_op", test_mat_mm_op },
   { "mat_linsolve", test_mat_linsolve },
   { "nmat_base", test_nmat_base },
+  { "nmat_new", test_nmat_new },
   { "nmat_linsolve", test_nmat_linsolve },
   { "nmat_linsolve_pa", test_nmat_linsolve_pa },
   { "mv", test_mat_mv},
diff --git a/gyrokinetic/unit/ctest_nodal_ops.c b/core/unit/ctest_nodal_ops.c
similarity index 50%
rename from gyrokinetic/unit/ctest_nodal_ops.c
rename to core/unit/ctest_nodal_ops.c
index fc248a31c0..6cf008ca6a 100644
--- a/gyrokinetic/unit/ctest_nodal_ops.c
+++ b/core/unit/ctest_nodal_ops.c
@@ -10,7 +10,6 @@
 #include <gkyl_util.h>
 #include <gkyl_basis.h>
 #include <gkyl_nodal_ops.h>
-#include <gkyl_deflate_zsurf.h>
 
 
 void check_same(struct gkyl_range range, struct gkyl_basis basis, struct gkyl_array *field1, struct gkyl_array* field2)
@@ -59,12 +58,12 @@ test_p1_2x(){
 #else
   struct gkyl_basis *basis_on_dev = &basis;
 #endif
-  
+
   // Project initial function
   struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   struct gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func, 0);
   gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
 #ifdef GKYL_HAVE_CUDA
   struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   gkyl_array_copy(funcdg_dev, funcdg);
@@ -103,7 +102,7 @@ test_p1_2x(){
   gkyl_array_copy(funcdg2, funcdg2_dev);
 #endif
   check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
 
 
 
@@ -142,12 +141,12 @@ test_p1_3x(){
 #else
   struct gkyl_basis *basis_on_dev = &basis;
 #endif
-  
+
   // Project initial function
   struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   struct gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func3d, 0);
   gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func3d.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func3d.gkyl");
 #ifdef GKYL_HAVE_CUDA
   struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   gkyl_array_copy(funcdg_dev, funcdg);
@@ -186,7 +185,7 @@ test_p1_3x(){
   gkyl_array_copy(funcdg2, funcdg2_dev);
 #endif
   check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func3d_2.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func3d_2.gkyl");
 
 
 
@@ -221,12 +220,12 @@ test_p1_interior_2x(){
   gkyl_cart_modal_serendip(&basis, 2, poly_order);
   bool use_gpu = false;
   struct gkyl_basis *basis_on_dev = &basis;
-  
+
   // Project initial function
   struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   struct gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func, 0);
   gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
   struct gkyl_array *funcdg_dev = funcdg;
 
   // Construct nrange and nodal field
@@ -249,7 +248,7 @@ test_p1_interior_2x(){
   gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, true);
 
   check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
 
 
 
@@ -276,12 +275,12 @@ test_p1_interior_3x(){
   gkyl_cart_modal_serendip(&basis, 3, poly_order);
   bool use_gpu = false;
   struct gkyl_basis *basis_on_dev = &basis;
-  
+
   // Project initial function
   struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   struct gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func3d, 0);
   gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func3d.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func3d.gkyl");
   struct gkyl_array *funcdg_dev = funcdg;
 
   // Construct nrange and nodal field
@@ -304,282 +303,18 @@ test_p1_interior_3x(){
   gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, true);
 
   check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func3d_2.gkyl");
-
-
-
-  gkyl_array_release(funcdg);
-  gkyl_array_release(funcdg2);
-  gkyl_array_release(nodal_fld);
-  gkyl_eval_on_nodes_release(eon);
-  gkyl_nodal_ops_release(n2m);
-
-}
-
-
-
-
-void
-test_p1_deflated(){
-  // create  grid, ranges, basis
-  double lower[] = { 0.0, -1.5 }, upper[] = { 1.5, 1.5 };
-  int cells[] = { 8, 16 };
-  struct gkyl_rect_grid grid;
-  gkyl_rect_grid_init(&grid, 2, lower, upper, cells);
-  struct gkyl_range local, local_ext;
-  int nghost[GKYL_MAX_CDIM] = { 1, 1 };
-  gkyl_create_grid_ranges(&grid, nghost, &local_ext, &local);
-  int poly_order = 1;
-  struct gkyl_basis basis;
-  gkyl_cart_modal_serendip(&basis, 2, poly_order);
-  bool use_gpu = false;
-#ifdef GKYL_HAVE_CUDA
-  use_gpu = true;
-  struct gkyl_basis *basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
-  gkyl_cart_modal_serendip_cu_dev(basis_on_dev, 2, poly_order);
-#else
-  struct gkyl_basis *basis_on_dev = &basis;
-#endif
-
-  // Project Initial Function
-  struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-  gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func, 0);
-  gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-  gkyl_array_copy(funcdg_dev, funcdg);
-#else
-  struct gkyl_array *funcdg_dev = funcdg;
-#endif
-
-  // Construct Nodal range, field, and nodal ops object
-  int nodes[3] = { 1, 1, 1 };
-  for (int d=0; d<grid.ndim; ++d)
-    nodes[d] = grid.cells[d] + 1;
-  struct gkyl_range nrange;
-  gkyl_range_init_from_shape(&nrange, grid.ndim, nodes);
-  struct gkyl_array* nodal_fld = gkyl_array_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *nodal_fld_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
-#else
-  struct gkyl_array *nodal_fld_dev = nodal_fld;
-#endif
-  struct gkyl_nodal_ops *n2m = gkyl_nodal_ops_new(&basis, &grid, use_gpu);
-
-  // Create deflated 1d grid, ranges, basis, nodal range
-  struct gkyl_rect_grid deflated_grid;
-  struct gkyl_basis deflated_basis;
-  struct gkyl_range deflated_local, deflated_local_ext;
-  double deflated_lower[1] = { grid.lower[0]}, deflated_upper[1] = { grid.upper[0]};
-  int deflated_cells[1] = { grid.cells[0] };
-  gkyl_rect_grid_init(&deflated_grid, 1, deflated_lower, deflated_upper, deflated_cells);
-  int deflated_nghost[GKYL_MAX_CDIM] = { 1 };
-  gkyl_create_grid_ranges(&deflated_grid, deflated_nghost, &deflated_local_ext, &deflated_local);
-  gkyl_cart_modal_serendip(&deflated_basis, 1, poly_order);
-#ifdef GKYL_HAVE_CUDA
-  use_gpu = true;
-  struct gkyl_basis *deflated_basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
-  gkyl_cart_modal_serendip_cu_dev(deflated_basis_on_dev, deflated_grid.ndim, poly_order);
-#else
-  struct gkyl_basis *deflated_basis_on_dev = &deflated_basis;
-#endif
-  int deflated_nodes[1];
-  for (int d=0; d<deflated_grid.ndim; ++d)
-    deflated_nodes[d] = deflated_grid.cells[d] + 1;
-  struct gkyl_range deflated_nrange;
-  gkyl_range_init_from_shape(&deflated_nrange, deflated_grid.ndim, deflated_nodes);
-
-  // Create deflated modal field 
-  struct gkyl_array *deflated_funcdg = gkyl_array_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *deflated_funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
-#else
-  struct gkyl_array *deflated_funcdg_dev = deflated_funcdg;
-#endif
-
-  // Created deflator objects and 1d nodal ops object
-  struct gkyl_deflate_zsurf *deflator_lo = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 0, use_gpu);
-  struct gkyl_deflate_zsurf *deflator_up = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 1, use_gpu);
-  struct gkyl_nodal_ops *n2m_1d = gkyl_nodal_ops_new(&deflated_basis, &deflated_grid, use_gpu);
-
-  // Loop along last dim (z) and call the 1d deflator to populate the 2d nodal field
-  for(int zidx = local.lower[1]; zidx <= local.upper[1]; zidx++){
-    gkyl_deflate_zsurf_advance(deflator_lo, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
-    gkyl_nodal_ops_m2n_deflated(n2m_1d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx-1);
-    if (zidx == local.upper[1]) {
-      gkyl_deflate_zsurf_advance(deflator_up, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
-      gkyl_nodal_ops_m2n_deflated(n2m_1d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx);
-    }
-  }
-
-  // Create output field
-  struct gkyl_array *funcdg2 = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *funcdg2_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-#else
-  struct gkyl_array *funcdg2_dev = funcdg2;
-#endif
-
-  // Transform back to modal
-  gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, false);
-#ifdef GKYL_HAVE_CUDA
-  gkyl_array_copy(funcdg2, funcdg2_dev);
-#endif
-  check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
-
-#ifdef GKYL_HAVE_CUDA
-  gkyl_cu_free(basis_on_dev);
-  gkyl_array_release(funcdg_dev);
-  gkyl_array_release(deflated_funcdg_dev);
-  gkyl_array_release(funcdg2_dev);
-  gkyl_array_release(nodal_fld_dev);
-#endif
-  gkyl_array_release(funcdg);
-  gkyl_array_release(deflated_funcdg);
-  gkyl_array_release(funcdg2);
-  gkyl_array_release(nodal_fld);
-
-  gkyl_eval_on_nodes_release(eon);
-  gkyl_nodal_ops_release(n2m);
-  gkyl_nodal_ops_release(n2m_1d);
-  gkyl_deflate_zsurf_release(deflator_lo);
-  gkyl_deflate_zsurf_release(deflator_up);
-
-}
-
-void test_p1_deflated_3d(){
-  // create  grid, ranges, basis
-  double lower[] = { 0.0, -1.5, -2 }, upper[] = { 1.5, 1.5,2 };
-  int cells[] = { 8, 16 , 12};
-  struct gkyl_rect_grid grid;
-  gkyl_rect_grid_init(&grid, 3, lower, upper, cells);
-  struct gkyl_range local, local_ext;
-  int nghost[GKYL_MAX_CDIM] = { 1, 1 ,1};
-  gkyl_create_grid_ranges(&grid, nghost, &local_ext, &local);
-  int poly_order = 1;
-  struct gkyl_basis basis;
-  gkyl_cart_modal_serendip(&basis, 3, poly_order);
-  bool use_gpu = false;
-#ifdef GKYL_HAVE_CUDA
-  use_gpu = true;
-  struct gkyl_basis *basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
-  gkyl_cart_modal_serendip_cu_dev(basis_on_dev, 3, poly_order);
-#else
-  struct gkyl_basis *basis_on_dev = &basis;
-#endif
-
-  // Project Initial Function
-  struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-  gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func3d, 0);
-  gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-  gkyl_array_copy(funcdg_dev, funcdg);
-#else
-  struct gkyl_array *funcdg_dev = funcdg;
-#endif
-
-  // Construct Nodal range, field, and nodal ops object
-  int nodes[3] = { 1, 1, 1 };
-  for (int d=0; d<grid.ndim; ++d)
-    nodes[d] = grid.cells[d] + 1;
-  struct gkyl_range nrange;
-  gkyl_range_init_from_shape(&nrange, grid.ndim, nodes);
-  struct gkyl_array* nodal_fld = gkyl_array_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *nodal_fld_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
-#else
-  struct gkyl_array *nodal_fld_dev = nodal_fld;
-#endif
-  struct gkyl_nodal_ops *n2m = gkyl_nodal_ops_new(&basis, &grid, use_gpu);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func3d_2.gkyl");
 
-  // Create deflated 2d grid, ranges, basis, nodal range
-  struct gkyl_rect_grid deflated_grid;
-  struct gkyl_basis deflated_basis;
-  struct gkyl_range deflated_local, deflated_local_ext;
-  double deflated_lower[2] = { grid.lower[0], grid.lower[1] };
-  double deflated_upper[2] = { grid.upper[0], grid.upper[1] };
-  int deflated_cells[2] = { grid.cells[0], grid.cells[1] };
-  gkyl_rect_grid_init(&deflated_grid, 2, deflated_lower, deflated_upper, deflated_cells);
-  int deflated_nghost[GKYL_MAX_CDIM] = { 1, 1 };
-  gkyl_create_grid_ranges(&deflated_grid, deflated_nghost, &deflated_local_ext, &deflated_local);
-  gkyl_cart_modal_serendip(&deflated_basis, 2, poly_order);
-#ifdef GKYL_HAVE_CUDA
-  use_gpu = true;
-  struct gkyl_basis *deflated_basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
-  gkyl_cart_modal_serendip_cu_dev(deflated_basis_on_dev, deflated_grid.ndim, poly_order);
-#else
-  struct gkyl_basis *deflated_basis_on_dev = &deflated_basis;
-#endif
-  int deflated_nodes[2];
-  for (int d=0; d<deflated_grid.ndim; ++d)
-    deflated_nodes[d] = deflated_grid.cells[d] + 1;
-  struct gkyl_range deflated_nrange;
-  gkyl_range_init_from_shape(&deflated_nrange, deflated_grid.ndim, deflated_nodes);
-
-  // Create deflated modal field 
-  struct gkyl_array *deflated_funcdg = gkyl_array_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *deflated_funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
-#else
-  struct gkyl_array *deflated_funcdg_dev = deflated_funcdg;
-#endif
 
-  // Created deflator objects and 1d nodal ops object
-  struct gkyl_deflate_zsurf *deflator_lo = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 0, use_gpu);
-  struct gkyl_deflate_zsurf *deflator_up = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 1, use_gpu);
-  struct gkyl_nodal_ops *n2m_2d = gkyl_nodal_ops_new(&deflated_basis, &deflated_grid, use_gpu);
-
-  // Loop along last dim (z) and call the 1d deflator to populate the 2d nodal field
-  for(int zidx = local.lower[2]; zidx <= local.upper[2]; zidx++){
-    gkyl_deflate_zsurf_advance(deflator_lo, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
-    gkyl_nodal_ops_m2n_deflated(n2m_2d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx-1);
-    if (zidx == local.upper[2]) {
-      gkyl_deflate_zsurf_advance(deflator_up, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
-      gkyl_nodal_ops_m2n_deflated(n2m_2d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx);
-    }
-  }
 
-  // Create output field
-  struct gkyl_array *funcdg2 = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-#ifdef GKYL_HAVE_CUDA
-  struct gkyl_array *funcdg2_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-#else
-  struct gkyl_array *funcdg2_dev = funcdg2;
-#endif
-
-  // Transform back to modal
-  gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, false);
-#ifdef GKYL_HAVE_CUDA
-  gkyl_array_copy(funcdg2, funcdg2_dev);
-#endif
-  check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
-
-#ifdef GKYL_HAVE_CUDA
-  gkyl_cu_free(basis_on_dev);
-  gkyl_array_release(funcdg_dev);
-  gkyl_array_release(deflated_funcdg_dev);
-  gkyl_array_release(funcdg2_dev);
-  gkyl_array_release(nodal_fld_dev);
-#endif
   gkyl_array_release(funcdg);
-  gkyl_array_release(deflated_funcdg);
   gkyl_array_release(funcdg2);
   gkyl_array_release(nodal_fld);
-
   gkyl_eval_on_nodes_release(eon);
   gkyl_nodal_ops_release(n2m);
-  gkyl_nodal_ops_release(n2m_2d);
-  gkyl_deflate_zsurf_release(deflator_lo);
-  gkyl_deflate_zsurf_release(deflator_up);
 
 }
 
-
 void
 test_p2_btype(enum gkyl_basis_type basis_type){
   // create  grid
@@ -606,14 +341,14 @@ test_p2_btype(enum gkyl_basis_type basis_type){
     gkyl_cart_modal_tensor(&basis, 2, poly_order);
 
   struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
-  
+
   gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func, 0);
   gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
   gkyl_eval_on_nodes_release(eon);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
 
   int nodes[3] = { 1, 1, 1 };
-                   
+
   for (int d=0; d<grid.ndim; ++d)
     nodes[d] = 2*(grid.cells[d]) + 1;
 
@@ -629,12 +364,12 @@ test_p2_btype(enum gkyl_basis_type basis_type){
   gkyl_nodal_ops_n2m(n2m, &basis, &grid, &nrange, &local, 1, nodal_fld, funcdg2, false);
   gkyl_nodal_ops_release(n2m);
   check_same(local, basis, funcdg, funcdg2);
-  gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
 
   gkyl_array_release(funcdg);
   gkyl_array_release(funcdg2);
   gkyl_array_release(nodal_fld);
-  
+
 }
 
 void
@@ -653,10 +388,8 @@ test_p2_tensor()
 TEST_LIST = {
   { "test_p1_interior_2x", test_p1_interior_2x},
   { "test_p1_interior_3x", test_p1_interior_3x},
-  {"test_p1_2x", test_p1_2x},
+  { "test_p1_2x", test_p1_2x},
   { "test_p1_3x", test_p1_3x},
-  { "test_p1_deflated", test_p1_deflated},
-  { "test_p1_deflated_3d", test_p1_deflated_3d},
   { "test_p2_ser", test_p2_ser},
   { "test_p2_tensor", test_p2_tensor},
   { NULL, NULL },
diff --git a/core/unit/ctest_range.c b/core/unit/ctest_range.c
index e36c6ef042..fc3f51fe3c 100644
--- a/core/unit/ctest_range.c
+++ b/core/unit/ctest_range.c
@@ -1337,6 +1337,128 @@ test_range_edge_match(void)
   TEST_CHECK( dir_ed.eloc == GKYL_NO_EDGE );
 }
 
+void
+test_range_init_shape()
+{
+  int lower[] = { 1, 1 }, upper[] = { 4, 5 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  TEST_CHECK( rng.ndim == 2 );
+  TEST_CHECK( rng.volume == 4*5 );
+  TEST_CHECK( gkyl_range_shape(&rng, 0) == 4 );
+  TEST_CHECK( gkyl_range_shape(&rng, 1) == 5 );
+  TEST_CHECK( !gkyl_range_is_sub_range(&rng) );
+}
+
+void
+test_range_init_from_shape_basic()
+{
+  int shape[] = { 3, 4, 5 };
+  struct gkyl_range rng;
+  gkyl_range_init_from_shape(&rng, 3, shape);
+
+  TEST_CHECK( rng.ndim == 3 );
+  TEST_CHECK( rng.volume == 60 );
+  TEST_CHECK( gkyl_range_shape(&rng, 0) == 3 );
+  TEST_CHECK( gkyl_range_shape(&rng, 1) == 4 );
+  TEST_CHECK( gkyl_range_shape(&rng, 2) == 5 );
+  // default lower index is 0
+  TEST_CHECK( rng.lower[0] == 0 );
+}
+
+void
+test_range_index_roundtrip()
+{
+  int lower[] = { 1, 1 }, upper[] = { 3, 4 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  // Every linear index in [0, volume) maps back to a valid multi-index
+  // whose forward index recovers the same linear index.
+  for (long loc=0; loc<rng.volume; ++loc) {
+    int idx[2];
+    gkyl_range_inv_idx(&rng, loc, idx);
+    TEST_CHECK( idx[0] >= rng.lower[0] && idx[0] <= rng.upper[0] );
+    TEST_CHECK( idx[1] >= rng.lower[1] && idx[1] <= rng.upper[1] );
+    long lin = gkyl_range_idx(&rng, idx);
+    TEST_CHECK( lin == loc );
+  }
+}
+
+void
+test_range_index_unique()
+{
+  int lower[] = { 0, 0 }, upper[] = { 2, 2 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  // All linear indices over the range are distinct and cover [0, volume).
+  int seen[9] = { 0 };
+  for (int i=lower[0]; i<=upper[0]; ++i) {
+    for (int j=lower[1]; j<=upper[1]; ++j) {
+      int idx[] = { i, j };
+      long lin = gkyl_range_idx(&rng, idx);
+      TEST_CHECK( lin >= 0 && lin < rng.volume );
+      seen[lin]++;
+    }
+  }
+  for (int k=0; k<9; ++k)
+    TEST_CHECK( seen[k] == 1 );
+}
+
+void
+test_range_iter_basic()
+{
+  int lower[] = { 1, 1 }, upper[] = { 3, 3 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &rng);
+
+  long count = 0;
+  while (gkyl_range_iter_next(&iter)) {
+    TEST_CHECK( iter.idx[0] >= 1 && iter.idx[0] <= 3 );
+    TEST_CHECK( iter.idx[1] >= 1 && iter.idx[1] <= 3 );
+    count++;
+  }
+  TEST_CHECK( count == rng.volume );
+  TEST_CHECK( count == 9 );
+}
+
+void
+test_range_shorten()
+{
+  int lower[] = { 1, 1 }, upper[] = { 4, 6 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  struct gkyl_range srng;
+  gkyl_range_shorten_from_above(&srng, &rng, 1, 1);
+  // Shortened to 1 cell in direction 1.
+  TEST_CHECK( gkyl_range_shape(&srng, 1) == 1 );
+  TEST_CHECK( gkyl_range_shape(&srng, 0) == 4 );
+  TEST_CHECK( srng.volume == 4 );
+}
+
+void
+test_sub_range_basic()
+{
+  int lower[] = { 1, 1 }, upper[] = { 10, 10 };
+  struct gkyl_range rng;
+  gkyl_range_init(&rng, 2, lower, upper);
+
+  int sublo[] = { 2, 2 }, subup[] = { 5, 5 };
+  struct gkyl_range sub;
+  gkyl_sub_range_init(&sub, &rng, sublo, subup);
+
+  TEST_CHECK( gkyl_range_is_sub_range(&sub) );
+  TEST_CHECK( sub.volume == 16 );
+  TEST_CHECK( gkyl_range_shape(&sub, 0) == 4 );
+  TEST_CHECK( gkyl_range_shape(&sub, 1) == 4 );
+}
+
 // CUDA specific tests
 #ifdef GKYL_HAVE_CUDA
 
@@ -1402,6 +1524,13 @@ TEST_LIST = {
   { "skin_ghost", test_skin_ghost },
   { "skin_ghost_with_corners", test_skin_ghost_with_corners },
   { "range_edge_match", test_range_edge_match },
+  { "range_init_shape", test_range_init_shape },
+  { "range_init_from_shape_basic", test_range_init_from_shape_basic },
+  { "range_index_roundtrip", test_range_index_roundtrip },
+  { "range_index_unique", test_range_index_unique },
+  { "range_iter_basic", test_range_iter_basic },
+  { "range_shorten", test_range_shorten },
+  { "sub_range_basic", test_sub_range_basic },
 #ifdef GKYL_HAVE_CUDA
   { "cu_range", test_cu_range },
 #endif  
diff --git a/core/unit/ctest_rect_grid.c b/core/unit/ctest_rect_grid.c
index de3412f551..735178be54 100644
--- a/core/unit/ctest_rect_grid.c
+++ b/core/unit/ctest_rect_grid.c
@@ -3,6 +3,7 @@
 #include <gkyl_alloc.h>
 #include <gkyl_rect_grid.h>
 #include <gkyl_util.h>
+#include <math.h>
 
 void test_grid_2d()
 {
@@ -215,6 +216,138 @@ void test_grid_io()
   TEST_CHECK( grid.cellVolume == grid2.cellVolume );
 }
 
+void
+test_grid_init_1d()
+{
+  double lower[] = { 0.0 }, upper[] = { 1.0 };
+  int cells[] = { 10 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+
+  TEST_CHECK( grid.ndim == 1 );
+  TEST_CHECK( grid.cells[0] == 10 );
+  TEST_CHECK( gkyl_compare_double(grid.lower[0], 0.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(grid.upper[0], 1.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(grid.dx[0], 0.1, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(grid.cellVolume, 0.1, 1e-15) );
+}
+
+void
+test_grid_init_3d()
+{
+  double lower[] = { -1.0, 0.0, 2.0 }, upper[] = { 1.0, 4.0, 6.0 };
+  int cells[] = { 4, 8, 2 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 3, lower, upper, cells);
+
+  TEST_CHECK( grid.ndim == 3 );
+  TEST_CHECK( gkyl_compare_double(grid.dx[0], 0.5, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(grid.dx[1], 0.5, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(grid.dx[2], 2.0, 1e-15) );
+  // cellVolume = product of dx
+  TEST_CHECK( gkyl_compare_double(grid.cellVolume, 0.5*0.5*2.0, 1e-15) );
+}
+
+void
+test_grid_cell_center()
+{
+  double lower[] = { 0.0, 0.0 }, upper[] = { 2.0, 2.0 };
+  int cells[] = { 2, 2 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 2, lower, upper, cells);
+
+  // dx = 1.0 in both directions; cell 1 center should be at 0.5
+  int idx[] = { 1, 1 };
+  double xc[2];
+  gkyl_rect_grid_cell_center(&grid, idx, xc);
+  TEST_CHECK( gkyl_compare_double(xc[0], 0.5, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(xc[1], 0.5, 1e-15) );
+
+  int idx2[] = { 2, 2 };
+  gkyl_rect_grid_cell_center(&grid, idx2, xc);
+  TEST_CHECK( gkyl_compare_double(xc[0], 1.5, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(xc[1], 1.5, 1e-15) );
+}
+
+void
+test_grid_ll_node()
+{
+  double lower[] = { 0.0 }, upper[] = { 1.0 };
+  int cells[] = { 4 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+
+  // dx = 0.25; lower-left node of cell i is lower + (i-1)*dx
+  int idx[] = { 1 };
+  double xc[1];
+  gkyl_rect_grid_ll_node(&grid, idx, xc);
+  TEST_CHECK( gkyl_compare_double(xc[0], 0.0, 1e-15) );
+
+  idx[0] = 3;
+  gkyl_rect_grid_ll_node(&grid, idx, xc);
+  TEST_CHECK( gkyl_compare_double(xc[0], 0.5, 1e-15) );
+
+  // cell center is ll_node + dx/2
+  double cc[1];
+  gkyl_rect_grid_cell_center(&grid, idx, cc);
+  TEST_CHECK( gkyl_compare_double(cc[0], xc[0] + 0.5*grid.dx[0], 1e-15) );
+}
+
+void
+test_grid_extents()
+{
+  double lower[] = { 0.0, 0.0 }, upper[] = { 1.0, 1.0 };
+  int cells[] = { 5, 7 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 2, lower, upper, cells);
+
+  int ext[2];
+  gkyl_rect_grid_extents(&grid, 0, ext);
+  TEST_CHECK( ext[0] == 1 );
+  TEST_CHECK( ext[1] == 5 );
+
+  gkyl_rect_grid_extents(&grid, 1, ext);
+  TEST_CHECK( ext[0] == 1 );
+  TEST_CHECK( ext[1] == 7 );
+}
+
+void
+test_grid_coord_idx()
+{
+  double lower[] = { 0.0 }, upper[] = { 10.0 };
+  int cells[] = { 10 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+
+  // dx = 1.0. x=0.5 -> cell 1, x=3.5 -> cell 4, x=9.9 -> cell 10
+  double xn[1];
+  int idx[1];
+
+  xn[0] = 0.5; gkyl_rect_grid_coord_idx(&grid, xn, idx);
+  TEST_CHECK( idx[0] == 1 );
+
+  xn[0] = 3.5; gkyl_rect_grid_coord_idx(&grid, xn, idx);
+  TEST_CHECK( idx[0] == 4 );
+
+  xn[0] = 9.9; gkyl_rect_grid_coord_idx(&grid, xn, idx);
+  TEST_CHECK( idx[0] == 10 );
+}
+
+void
+test_grid_cmp()
+{
+  double lower[] = { 0.0, 0.0 }, upper[] = { 1.0, 1.0 };
+  int cells[] = { 4, 4 };
+  struct gkyl_rect_grid g1, g2, g3;
+  gkyl_rect_grid_init(&g1, 2, lower, upper, cells);
+  gkyl_rect_grid_init(&g2, 2, lower, upper, cells);
+  TEST_CHECK( gkyl_rect_grid_cmp(&g1, &g2) );
+
+  int cells2[] = { 4, 8 };
+  gkyl_rect_grid_init(&g3, 2, lower, upper, cells2);
+  TEST_CHECK( !gkyl_rect_grid_cmp(&g1, &g3) );
+}
+
 // CUDA specific tests
 #ifdef GKYL_HAVE_CUDA
 
@@ -239,6 +372,13 @@ TEST_LIST = {
   { "grid_find_cell_2d", test_find_cell_2d },
   { "grid_find_cell_3d", test_find_cell_3d },
   { "grid_io", test_grid_io },
+  { "grid_init_1d", test_grid_init_1d },
+  { "grid_init_3d", test_grid_init_3d },
+  { "grid_cell_center", test_grid_cell_center },
+  { "grid_ll_node", test_grid_ll_node },
+  { "grid_extents", test_grid_extents },
+  { "grid_coord_idx", test_grid_coord_idx },
+  { "grid_cmp", test_grid_cmp },
 #ifdef GKYL_HAVE_CUDA
   { "cu_grid_2d", test_cu_grid_2d },
 #endif  
diff --git a/core/unit/ctest_util.c b/core/unit/ctest_util.c
new file mode 100644
index 0000000000..105f5af2ba
--- /dev/null
+++ b/core/unit/ctest_util.c
@@ -0,0 +1,186 @@
+// Unit tests for utility helpers in gkyl_util.h / util.c and constants in gkyl_const.h
+#include <acutest.h>
+#include <gkyl_util.h>
+#include <gkyl_const.h>
+#include <math.h>
+
+void
+test_compare_double()
+{
+  TEST_CHECK( gkyl_compare_double(1.0, 1.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(0.0, 0.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(1.0, 1.0 + 1e-16, 1e-12) );
+  TEST_CHECK( !gkyl_compare_double(1.0, 2.0, 1e-12) );
+  TEST_CHECK( !gkyl_compare_double(1.0, 1.1, 1e-6) );
+  // NaN never compares equal
+  TEST_CHECK( !gkyl_compare_double(NAN, 1.0, 1e-6) );
+  TEST_CHECK( !gkyl_compare_double(1.0, NAN, 1e-6) );
+  // near zero
+  TEST_CHECK( gkyl_compare_double(0.0, 1e-20, 1e-12) );
+  // symmetric
+  TEST_CHECK( gkyl_compare_double(3.0, 3.0+1e-13, 1e-9) == gkyl_compare_double(3.0+1e-13, 3.0, 1e-9) );
+}
+
+void
+test_compare_float()
+{
+  TEST_CHECK( gkyl_compare_float(1.0f, 1.0f, 1e-6f) );
+  TEST_CHECK( !gkyl_compare_float(1.0f, 2.0f, 1e-6f) );
+  TEST_CHECK( gkyl_compare_float(0.0f, 0.0f, 1e-6f) );
+  TEST_CHECK( gkyl_compare_float(100.0f, 100.0f, 1e-6f) );
+}
+
+void
+test_copy_int_arr()
+{
+  int in[5] = { 1, 2, 3, 4, 5 };
+  int out[5] = { 0 };
+  gkyl_copy_int_arr(5, in, out);
+  for (int i=0; i<5; ++i)
+    TEST_CHECK( out[i] == in[i] );
+
+  // partial copy
+  int out2[5] = { 9, 9, 9, 9, 9 };
+  gkyl_copy_int_arr(3, in, out2);
+  TEST_CHECK( out2[0]==1 && out2[1]==2 && out2[2]==3 );
+  TEST_CHECK( out2[3]==9 && out2[4]==9 );
+}
+
+void
+test_copy_long_arr()
+{
+  long in[4] = { 10L, 20L, 30L, 40L };
+  long out[4] = { 0 };
+  gkyl_copy_long_arr(4, in, out);
+  for (int i=0; i<4; ++i)
+    TEST_CHECK( out[i] == in[i] );
+}
+
+void
+test_copy_double_arr()
+{
+  double in[4] = { 1.5, -2.5, 3.25, 0.0 };
+  double out[4] = { 0 };
+  gkyl_copy_double_arr(4, in, out);
+  for (int i=0; i<4; ++i)
+    TEST_CHECK( out[i] == in[i] );
+}
+
+void
+test_int_div_up()
+{
+  TEST_CHECK( gkyl_int_div_up(10, 5) == 2 );
+  TEST_CHECK( gkyl_int_div_up(11, 5) == 3 );
+  TEST_CHECK( gkyl_int_div_up(9, 5) == 2 );
+  TEST_CHECK( gkyl_int_div_up(0, 5) == 0 );
+  TEST_CHECK( gkyl_int_div_up(5, 5) == 1 );
+  TEST_CHECK( gkyl_int_div_up(1, 5) == 1 );
+  TEST_CHECK( gkyl_int_div_up(100, 7) == 15 );
+}
+
+void
+test_minmod_util()
+{
+  // all same sign positive -> min
+  TEST_CHECK( gkyl_minmod(2.0, 3.0, 4.0) == 2.0 );
+  // all negative -> max (closest to zero)
+  TEST_CHECK( gkyl_minmod(-2.0, -3.0, -4.0) == -2.0 );
+  // mixed signs -> 0
+  TEST_CHECK( gkyl_minmod(2.0, -3.0, 4.0) == 0.0 );
+  TEST_CHECK( gkyl_minmod(-1.0, 2.0, 3.0) == 0.0 );
+}
+
+void
+test_sgn_macro()
+{
+  TEST_CHECK( GKYL_SGN(5.0) == 1.0 );
+  TEST_CHECK( GKYL_SGN(-5.0) == -1.0 );
+  TEST_CHECK( GKYL_SGN(0.0) == 1.0 );
+}
+
+void
+test_minmax_macros()
+{
+  TEST_CHECK( GKYL_MIN2(3, 5) == 3 );
+  TEST_CHECK( GKYL_MIN2(5, 3) == 3 );
+  TEST_CHECK( GKYL_MAX2(3, 5) == 5 );
+  TEST_CHECK( GKYL_MAX2(5, 3) == 5 );
+  TEST_CHECK( GKYL_MIN2(-2.5, 1.0) == -2.5 );
+  TEST_CHECK( GKYL_MAX2(-2.5, 1.0) == 1.0 );
+}
+
+void
+test_tm_trigger()
+{
+  // Trigger every 1.0 time units.
+  struct gkyl_tm_trigger tmt = { .dt = 1.0, .tcurr = 0.0, .curr = 0 };
+
+  TEST_CHECK( gkyl_tm_trigger_check_and_bump(&tmt, 0.0) == 1 );
+  TEST_CHECK( tmt.curr == 1 );
+
+  // Below next threshold -> no trigger.
+  TEST_CHECK( gkyl_tm_trigger_check_and_bump(&tmt, 0.5) == 0 );
+  TEST_CHECK( tmt.curr == 1 );
+
+  TEST_CHECK( gkyl_tm_trigger_check_and_bump(&tmt, 1.0) == 1 );
+  TEST_CHECK( tmt.curr == 2 );
+
+  TEST_CHECK( gkyl_tm_trigger_check_and_bump(&tmt, 1.9) == 0 );
+  TEST_CHECK( gkyl_tm_trigger_check_and_bump(&tmt, 2.3) == 1 );
+  TEST_CHECK( tmt.curr == 3 );
+}
+
+void
+test_search_str_int_pair()
+{
+  struct gkyl_str_int_pair pairs[] = {
+    { "one", 1 },
+    { "two", 2 },
+    { "three", 3 },
+    { 0, 0 }
+  };
+
+  TEST_CHECK( gkyl_search_str_int_pair_by_str(pairs, "two", -1) == 2 );
+  TEST_CHECK( gkyl_search_str_int_pair_by_str(pairs, "one", -1) == 1 );
+  TEST_CHECK( gkyl_search_str_int_pair_by_str(pairs, "missing", -1) == -1 );
+
+  const char *s = gkyl_search_str_int_pair_by_int(pairs, 3, "none");
+  TEST_CHECK( strcmp(s, "three") == 0 );
+  const char *s2 = gkyl_search_str_int_pair_by_int(pairs, 99, "none");
+  TEST_CHECK( strcmp(s2, "none") == 0 );
+}
+
+void
+test_constants()
+{
+  // c = 1/sqrt(mu0 * eps0)
+  double c = 1.0/sqrt(GKYL_MU0*GKYL_EPSILON0);
+  TEST_CHECK( gkyl_compare_double(c, GKYL_SPEED_OF_LIGHT, 1e-6) );
+
+  // eV->Kelvin conversion factor consistency
+  double ev2k = GKYL_ELEMENTARY_CHARGE/GKYL_BOLTZMANN_CONSTANT;
+  TEST_CHECK( gkyl_compare_double(ev2k, GKYL_EV2KELVIN, 1e-9) );
+
+  // proton mass much larger than electron mass
+  TEST_CHECK( GKYL_PROTON_MASS > 1000.0*GKYL_ELECTRON_MASS );
+
+  // pi and e values
+  TEST_CHECK( gkyl_compare_double(GKYL_PI, M_PI, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(GKYL_E, M_E, 1e-14) );
+}
+
+TEST_LIST = {
+  { "compare_double", test_compare_double },
+  { "compare_float", test_compare_float },
+  { "copy_int_arr", test_copy_int_arr },
+  { "copy_long_arr", test_copy_long_arr },
+  { "copy_double_arr", test_copy_double_arr },
+  { "int_div_up", test_int_div_up },
+  { "minmod_util", test_minmod_util },
+  { "sgn_macro", test_sgn_macro },
+  { "minmax_macros", test_minmax_macros },
+  { "tm_trigger", test_tm_trigger },
+  { "search_str_int_pair", test_search_str_int_pair },
+  { "constants", test_constants },
+  { NULL, NULL },
+};
diff --git a/core/unit/ctest_vec3.c b/core/unit/ctest_vec3.c
new file mode 100644
index 0000000000..74d4fdcb72
--- /dev/null
+++ b/core/unit/ctest_vec3.c
@@ -0,0 +1,215 @@
+// Unit tests for the inline vector and limiter helpers in gkyl_math.h
+#include <acutest.h>
+#include <gkyl_math.h>
+#include <gkyl_util.h>
+#include <math.h>
+
+void
+test_vec3_basic()
+{
+  struct gkyl_vec3 z = gkyl_vec3_zeros();
+  TEST_CHECK( z.x[0] == 0.0 );
+  TEST_CHECK( z.x[1] == 0.0 );
+  TEST_CHECK( z.x[2] == 0.0 );
+
+  struct gkyl_vec3 a = gkyl_vec3_new(1.0, 2.0, 3.0);
+  TEST_CHECK( a.x[0] == 1.0 );
+  TEST_CHECK( a.x[1] == 2.0 );
+  TEST_CHECK( a.x[2] == 3.0 );
+}
+
+void
+test_vec3_scale()
+{
+  struct gkyl_vec3 a = gkyl_vec3_new(1.0, -2.0, 3.5);
+  struct gkyl_vec3 b = gkyl_vec3_scale(2.0, a);
+  TEST_CHECK( b.x[0] == 2.0 );
+  TEST_CHECK( b.x[1] == -4.0 );
+  TEST_CHECK( b.x[2] == 7.0 );
+
+  struct gkyl_vec3 c = gkyl_vec3_scale(0.0, a);
+  TEST_CHECK( c.x[0] == 0.0 && c.x[1] == 0.0 && c.x[2] == 0.0 );
+}
+
+void
+test_vec3_add_sub()
+{
+  struct gkyl_vec3 a = gkyl_vec3_new(1.0, 2.0, 3.0);
+  struct gkyl_vec3 b = gkyl_vec3_new(4.0, -1.0, 0.5);
+
+  struct gkyl_vec3 s = gkyl_vec3_add(a, b);
+  TEST_CHECK( s.x[0] == 5.0 );
+  TEST_CHECK( s.x[1] == 1.0 );
+  TEST_CHECK( s.x[2] == 3.5 );
+
+  struct gkyl_vec3 d = gkyl_vec3_sub(a, b);
+  TEST_CHECK( d.x[0] == -3.0 );
+  TEST_CHECK( d.x[1] == 3.0 );
+  TEST_CHECK( d.x[2] == 2.5 );
+
+  // a - a = 0
+  struct gkyl_vec3 z = gkyl_vec3_sub(a, a);
+  TEST_CHECK( z.x[0] == 0.0 && z.x[1] == 0.0 && z.x[2] == 0.0 );
+}
+
+void
+test_vec3_len_norm()
+{
+  struct gkyl_vec3 a = gkyl_vec3_new(3.0, 4.0, 0.0);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_len(a), 5.0, 1e-15) );
+
+  struct gkyl_vec3 n = gkyl_vec3_norm(a);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_len(n), 1.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(n.x[0], 0.6, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(n.x[1], 0.8, 1e-15) );
+
+  struct gkyl_vec3 b = gkyl_vec3_new(1.0, 2.0, 2.0);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_len(b), 3.0, 1e-15) );
+}
+
+void
+test_vec3_dot()
+{
+  struct gkyl_vec3 a = gkyl_vec3_new(1.0, 2.0, 3.0);
+  struct gkyl_vec3 b = gkyl_vec3_new(4.0, 5.0, 6.0);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_dot(a, b), 32.0, 1e-15) );
+
+  // orthogonal vectors
+  struct gkyl_vec3 e1 = gkyl_vec3_new(1.0, 0.0, 0.0);
+  struct gkyl_vec3 e2 = gkyl_vec3_new(0.0, 1.0, 0.0);
+  TEST_CHECK( gkyl_vec3_dot(e1, e2) == 0.0 );
+
+  // dot with self == len^2
+  double l = gkyl_vec3_len(a);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_dot(a, a), l*l, 1e-13) );
+}
+
+void
+test_vec3_cross()
+{
+  struct gkyl_vec3 e1 = gkyl_vec3_new(1.0, 0.0, 0.0);
+  struct gkyl_vec3 e2 = gkyl_vec3_new(0.0, 1.0, 0.0);
+  struct gkyl_vec3 e3 = gkyl_vec3_cross(e1, e2);
+  TEST_CHECK( gkyl_compare_double(e3.x[0], 0.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(e3.x[1], 0.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(e3.x[2], 1.0, 1e-15) );
+
+  // anti-commutativity: a x b = -(b x a)
+  struct gkyl_vec3 a = gkyl_vec3_new(2.0, -1.0, 3.0);
+  struct gkyl_vec3 b = gkyl_vec3_new(0.5, 4.0, -2.0);
+  struct gkyl_vec3 ab = gkyl_vec3_cross(a, b);
+  struct gkyl_vec3 ba = gkyl_vec3_cross(b, a);
+  TEST_CHECK( gkyl_compare_double(ab.x[0], -ba.x[0], 1e-14) );
+  TEST_CHECK( gkyl_compare_double(ab.x[1], -ba.x[1], 1e-14) );
+  TEST_CHECK( gkyl_compare_double(ab.x[2], -ba.x[2], 1e-14) );
+
+  // cross product is orthogonal to both inputs
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_dot(ab, a), 0.0, 1e-13) );
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_dot(ab, b), 0.0, 1e-13) );
+
+  // a x a = 0
+  struct gkyl_vec3 z = gkyl_vec3_cross(a, a);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_len(z), 0.0, 1e-14) );
+}
+
+void
+test_vec3_triple()
+{
+  struct gkyl_vec3 e1 = gkyl_vec3_new(1.0, 0.0, 0.0);
+  struct gkyl_vec3 e2 = gkyl_vec3_new(0.0, 1.0, 0.0);
+  struct gkyl_vec3 e3 = gkyl_vec3_new(0.0, 0.0, 1.0);
+  // triple product of unit basis is the determinant == 1
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_triple(e1, e2, e3), 1.0, 1e-15) );
+
+  // general: a.(b x c) equals scalar triple via determinant
+  struct gkyl_vec3 a = gkyl_vec3_new(1.0, 2.0, 3.0);
+  struct gkyl_vec3 b = gkyl_vec3_new(4.0, 5.0, 6.0);
+  struct gkyl_vec3 c = gkyl_vec3_new(7.0, 8.0, 10.0);
+  double det = 1.0*(5.0*10.0-6.0*8.0) - 2.0*(4.0*10.0-6.0*7.0) + 3.0*(4.0*8.0-5.0*7.0);
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_triple(a, b, c), det, 1e-12) );
+
+  // coplanar vectors give zero triple product
+  TEST_CHECK( gkyl_compare_double(gkyl_vec3_triple(e1, e2, gkyl_vec3_add(e1,e2)), 0.0, 1e-14) );
+}
+
+void
+test_vec3_polar()
+{
+  // At phi=0 the contravariant->cartesian transform should be near-identity in x.
+  struct gkyl_vec3 pin = gkyl_vec3_new(2.0, 0.0, 5.0);
+  struct gkyl_vec3 out = gkyl_vec3_polar_con_to_cart(1.0, 0.0, pin);
+  TEST_CHECK( gkyl_compare_double(out.x[0], 2.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(out.x[1], 0.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(out.x[2], 5.0, 1e-14) );
+
+  // con_to_cov scales the angular component by r^2 and leaves others.
+  struct gkyl_vec3 cov = gkyl_vec3_polar_con_to_cov(2.0, gkyl_vec3_new(1.0, 3.0, 7.0));
+  TEST_CHECK( gkyl_compare_double(cov.x[0], 1.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(cov.x[1], 3.0*4.0, 1e-14) );
+  TEST_CHECK( gkyl_compare_double(cov.x[2], 7.0, 1e-14) );
+}
+
+void
+test_minmod_2()
+{
+  TEST_CHECK( gkyl_minmod_2(2.0, 3.0) == 2.0 );
+  TEST_CHECK( gkyl_minmod_2(3.0, 2.0) == 2.0 );
+  TEST_CHECK( gkyl_minmod_2(-2.0, -3.0) == -2.0 );
+  TEST_CHECK( gkyl_minmod_2(-3.0, -2.0) == -2.0 );
+  // opposite signs -> 0
+  TEST_CHECK( gkyl_minmod_2(2.0, -3.0) == 0.0 );
+  TEST_CHECK( gkyl_minmod_2(-2.0, 3.0) == 0.0 );
+  TEST_CHECK( gkyl_minmod_2(0.0, 5.0) == 0.0 );
+}
+
+void
+test_minmod_3()
+{
+  TEST_CHECK( gkyl_minmod_3(2.0, 3.0, 4.0) == 2.0 );
+  TEST_CHECK( gkyl_minmod_3(-2.0, -3.0, -4.0) == -2.0 );
+  TEST_CHECK( gkyl_minmod_3(2.0, -3.0, 4.0) == 0.0 );
+  TEST_CHECK( gkyl_minmod_3(1.0, 2.0, 0.0) == 0.0 );
+}
+
+void
+test_minmod_4()
+{
+  TEST_CHECK( gkyl_minmod_4(2.0, 3.0, 4.0, 5.0) == 2.0 );
+  TEST_CHECK( gkyl_minmod_4(-2.0, -3.0, -4.0, -5.0) == -2.0 );
+  TEST_CHECK( gkyl_minmod_4(2.0, 3.0, -4.0, 5.0) == 0.0 );
+}
+
+void
+test_median()
+{
+  TEST_CHECK( gkyl_compare_double(gkyl_median(1.0, 2.0, 3.0), 2.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(gkyl_median(3.0, 1.0, 2.0), 2.0, 1e-15) );
+  TEST_CHECK( gkyl_compare_double(gkyl_median(2.0, 3.0, 1.0), 2.0, 1e-15) );
+}
+
+void
+test_min_max_3()
+{
+  TEST_CHECK( gkyl_min_3(1.0, 2.0, 3.0) == 1.0 );
+  TEST_CHECK( gkyl_min_3(3.0, -1.0, 2.0) == -1.0 );
+  TEST_CHECK( gkyl_max_3(1.0, 2.0, 3.0) == 3.0 );
+  TEST_CHECK( gkyl_max_3(3.0, -1.0, 2.0) == 3.0 );
+  TEST_CHECK( gkyl_max_3(-5.0, -2.0, -9.0) == -2.0 );
+}
+
+TEST_LIST = {
+  { "vec3_basic", test_vec3_basic },
+  { "vec3_scale", test_vec3_scale },
+  { "vec3_add_sub", test_vec3_add_sub },
+  { "vec3_len_norm", test_vec3_len_norm },
+  { "vec3_dot", test_vec3_dot },
+  { "vec3_cross", test_vec3_cross },
+  { "vec3_triple", test_vec3_triple },
+  { "vec3_polar", test_vec3_polar },
+  { "minmod_2", test_minmod_2 },
+  { "minmod_3", test_minmod_3 },
+  { "minmod_4", test_minmod_4 },
+  { "median", test_median },
+  { "min_max_3", test_min_max_3 },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_bc_basic_gyrokinetic.c b/gyrokinetic/unit/ctest_bc_basic_gyrokinetic.c
new file mode 100644
index 0000000000..08284cc982
--- /dev/null
+++ b/gyrokinetic/unit/ctest_bc_basic_gyrokinetic.c
@@ -0,0 +1,138 @@
+// Tests for the gyrokinetic basic-BC updater (gkyl_bc_basic_gyrokinetic).
+//
+// 1) Constructor field checks: dir/edge/cdim/bctype recorded, skin/ghost
+//    range pointers stored.
+// 2) COPY-BC advance compute check: fill the skin cell with known values,
+//    run the updater, and verify the ghost cell ends up holding an exact
+//    copy of the skin data (the COPY boundary condition).
+#include <acutest.h>
+
+#include <gkyl_array.h>
+#include <gkyl_array_ops.h>
+#include <gkyl_basis.h>
+#include <gkyl_bc_basic_gyrokinetic.h>
+#include <gkyl_bc_basic_gyrokinetic_priv.h>
+#include <gkyl_gk_bc_type.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_decomp.h>
+#include <gkyl_rect_grid.h>
+
+// Build a 1D conf grid + skin/ghost ranges at the given edge.
+static void
+setup_1x(enum gkyl_edge_loc edge, struct gkyl_basis *basis,
+  struct gkyl_range *local, struct gkyl_range *local_ext,
+  struct gkyl_range *skin_r, struct gkyl_range *ghost_r)
+{
+  int cells[] = { 8 };
+  int ghost[] = { 1 };
+  double lower[] = { 0.0 }, upper[] = { 1.0 };
+
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+  gkyl_create_grid_ranges(&grid, ghost, local_ext, local);
+
+  gkyl_cart_modal_serendip(basis, 1, 1);
+
+  gkyl_skin_ghost_ranges(skin_r, ghost_r, 0, edge, local_ext, ghost);
+}
+
+void
+test_bc_basic_ctor()
+{
+  struct gkyl_basis basis;
+  struct gkyl_range local, local_ext, skin_r, ghost_r;
+  setup_1x(GKYL_LOWER_EDGE, &basis, &local, &local_ext, &skin_r, &ghost_r);
+
+  struct gkyl_bc_basic_gyrokinetic *bc = gkyl_bc_basic_gyrokinetic_new(
+    0, GKYL_LOWER_EDGE, GKYL_BC_GK_SPECIES_COPY, &basis, &skin_r, &ghost_r,
+    basis.num_basis, 1, false);
+
+  TEST_CHECK( bc != NULL );
+  TEST_CHECK( bc->dir == 0 );
+  TEST_CHECK( bc->cdim == 1 );
+  TEST_CHECK( bc->edge == GKYL_LOWER_EDGE );
+  TEST_CHECK( bc->bctype == GKYL_BC_GK_SPECIES_COPY );
+  TEST_CHECK( bc->skin_r == &skin_r );
+  TEST_CHECK( bc->ghost_r == &ghost_r );
+  TEST_CHECK( bc->use_gpu == false );
+  TEST_CHECK( bc->array_copy_func != NULL );
+
+  gkyl_bc_basic_gyrokinetic_release(bc);
+}
+
+// Verify that applying a COPY BC fills the ghost cell with an exact copy of
+// the skin cell.
+static void
+check_copy_advance(enum gkyl_edge_loc edge)
+{
+  struct gkyl_basis basis;
+  struct gkyl_range local, local_ext, skin_r, ghost_r;
+  setup_1x(edge, &basis, &local, &local_ext, &skin_r, &ghost_r);
+
+  int nc = basis.num_basis; // 2 for 1x p1
+
+  struct gkyl_array *f = gkyl_array_new(GKYL_DOUBLE, nc, local_ext.volume);
+  gkyl_array_clear(f, 0.0);
+
+  // Buffer big enough to hold the skin cells.
+  long buff_sz = skin_r.volume;
+  struct gkyl_array *buff = gkyl_array_new(GKYL_DOUBLE, nc, buff_sz);
+  gkyl_array_clear(buff, 0.0);
+
+  // Fill the skin cell with known, distinct values.
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &skin_r);
+  while (gkyl_range_iter_next(&iter)) {
+    long sidx = gkyl_range_idx(&skin_r, iter.idx);
+    double *fs = gkyl_array_fetch(f, sidx);
+    for (int c=0; c<nc; c++)
+      fs[c] = 10.0 + c + 0.5;
+  }
+
+  struct gkyl_bc_basic_gyrokinetic *bc = gkyl_bc_basic_gyrokinetic_new(
+    0, edge, GKYL_BC_GK_SPECIES_COPY, &basis, &skin_r, &ghost_r, nc, 1, false);
+
+  gkyl_bc_basic_gyrokinetic_advance(bc, buff, f);
+
+  // Ghost cell should now equal the skin values.
+  gkyl_range_iter_init(&iter, &ghost_r);
+  while (gkyl_range_iter_next(&iter)) {
+    long gidx = gkyl_range_idx(&ghost_r, iter.idx);
+    const double *fg = gkyl_array_cfetch(f, gidx);
+    for (int c=0; c<nc; c++)
+      TEST_CHECK( gkyl_compare(fg[c], 10.0 + c + 0.5, 1e-14) );
+  }
+
+  gkyl_bc_basic_gyrokinetic_release(bc);
+  gkyl_array_release(f);
+  gkyl_array_release(buff);
+}
+
+void test_bc_basic_copy_lower() { check_copy_advance(GKYL_LOWER_EDGE); }
+void test_bc_basic_copy_upper() { check_copy_advance(GKYL_UPPER_EDGE); }
+
+void
+test_bc_basic_ctor_upper_reflect()
+{
+  struct gkyl_basis basis;
+  struct gkyl_range local, local_ext, skin_r, ghost_r;
+  setup_1x(GKYL_UPPER_EDGE, &basis, &local, &local_ext, &skin_r, &ghost_r);
+
+  struct gkyl_bc_basic_gyrokinetic *bc = gkyl_bc_basic_gyrokinetic_new(
+    0, GKYL_UPPER_EDGE, GKYL_BC_GK_SPECIES_REFLECT, &basis, &skin_r, &ghost_r,
+    basis.num_basis, 1, false);
+
+  TEST_CHECK( bc != NULL );
+  TEST_CHECK( bc->edge == GKYL_UPPER_EDGE );
+  TEST_CHECK( bc->bctype == GKYL_BC_GK_SPECIES_REFLECT );
+
+  gkyl_bc_basic_gyrokinetic_release(bc);
+}
+
+TEST_LIST = {
+  { "bc_basic_ctor", test_bc_basic_ctor },
+  { "bc_basic_copy_lower", test_bc_basic_copy_lower },
+  { "bc_basic_copy_upper", test_bc_basic_copy_upper },
+  { "bc_basic_ctor_upper_reflect", test_bc_basic_ctor_upper_reflect },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_bc_twistshift.c b/gyrokinetic/unit/ctest_bc_twistshift.c
index 707d48c149..5fd12e713d 100644
--- a/gyrokinetic/unit/ctest_bc_twistshift.c
+++ b/gyrokinetic/unit/ctest_bc_twistshift.c
@@ -358,7 +358,7 @@ test_bc_twistshift_3x_fig6_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig6_tar.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig6_tar.gkyl");
   }
 
   if (check_distf) {
@@ -409,7 +409,7 @@ test_bc_twistshift_3x_fig6_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig6_tar_shifted.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig6_tar_shifted.gkyl");
   }
 
   if (check_distf) {
@@ -615,7 +615,7 @@ test_bc_twistshift_3x2v_fig6_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig6_tar.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig6_tar.gkyl");
   }
 
   // Compute the integrated moments of the skin cell and the ghost cell.
@@ -750,7 +750,7 @@ test_bc_twistshift_3x2v_fig6_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig6_tar_shifted.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig6_tar_shifted.gkyl");
   }
 
   gkyl_dg_updater_moment_gyrokinetic_advance(mcalc,
@@ -1025,7 +1025,7 @@ test_bc_twistshift_3x_fig11_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig11_tar.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x_fig11_tar.gkyl");
   }
 
   if (check_distf) {
@@ -1780,8 +1780,8 @@ test_bc_twistshift_3x2v_fig11_wcells(const int *cells, enum gkyl_edge_loc edge,
       .basis_type = basis.id
     }
   );
-  if (write_f)
-    gkyl_grid_sub_array_write(&grid, &local, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig11_do.gkyl");
+  // if (write_f)
+  //   gkyl_grid_sub_array_write(&grid, &local, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig11_do.gkyl");
 
   // Create a range only extended in bc_dir.
   struct gkyl_range update_rng;
@@ -1842,7 +1842,7 @@ test_bc_twistshift_3x2v_fig11_wcells(const int *cells, enum gkyl_edge_loc edge,
     }
     struct gkyl_rect_grid grid_ext;
     gkyl_rect_grid_init(&grid_ext, ndim, lower_ext, upper_ext, cells_ext);
-    gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig11_tar.gkyl");
+    // gkyl_grid_sub_array_write(&grid_ext, &local_ext, mt, distf_ho, "ctest_bc_twistshift_3x2v_fig11_tar.gkyl");
   }
 
   // Compute the integrated moments of the skin cell and the ghost cell.
diff --git a/gyrokinetic/unit/ctest_block_tensor.c b/gyrokinetic/unit/ctest_block_tensor.c
index a1dd2891d2..bfb4ae8026 100644
--- a/gyrokinetic/unit/ctest_block_tensor.c
+++ b/gyrokinetic/unit/ctest_block_tensor.c
@@ -75,8 +75,8 @@ void test_cartesian_2x_onecell()
   gkyl_eval_on_nodes_advance(proj, 0.0, &local, dxdz);
   gkyl_eval_on_nodes_advance(proj, 0.0, &local, dzdx);
   gkyl_eval_on_nodes_release(proj);
-  gkyl_grid_sub_array_write(&grid, &local, 0, dxdz, "dxdz.gkyl");
-  gkyl_grid_sub_array_write(&grid, &local, 0, dzdx, "dzdx.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, dxdz, "dxdz.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, dzdx, "dzdx.gkyl");
 
   struct bc_block_tensor *bt = gkyl_bc_block_tensor_new(&grid, &local, &local_ext, &basis, false);
 
@@ -140,14 +140,14 @@ void test_cartesian_2x_z()
   gkyl_eval_on_nodes *proj2 = gkyl_eval_on_nodes_new(&grid2, &basis, 9, &proj_one, 0);
   gkyl_eval_on_nodes_advance(proj2, 0.0, &local2, dzdx2);
   gkyl_eval_on_nodes_release(proj2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
 
   // Block 1 tangents
   struct gkyl_array *dxdz1 = gkyl_array_new(GKYL_DOUBLE, 9*basis.num_basis, local_ext1.volume);
   gkyl_eval_on_nodes *proj1 = gkyl_eval_on_nodes_new(&grid1, &basis, 9, &proj_one, 0);
   gkyl_eval_on_nodes_advance(proj1, 0.0, &local1, dxdz1);
   gkyl_eval_on_nodes_release(proj1);
-  gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
+  // gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
 
 
 
@@ -156,7 +156,7 @@ void test_cartesian_2x_z()
   int edge2 = 0; //lower edge
   int dir = 1; // second direction
   gkyl_bc_block_tensor_advance(bt, dir, edge1, edge2, dxdz1, dzdx2, &local1, &local2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0,bt->tensor, "tji.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0,bt->tensor, "tji.gkyl");
 
   gkyl_array_release(dxdz1);
   gkyl_array_release(dzdx2);
@@ -195,14 +195,14 @@ void test_cartesian_2x_x()
   gkyl_eval_on_nodes *proj2 = gkyl_eval_on_nodes_new(&grid2, &basis, 9, &proj_one, 0);
   gkyl_eval_on_nodes_advance(proj2, 0.0, &local2, dzdx2);
   gkyl_eval_on_nodes_release(proj2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
 
   // Block 1 tangents
   struct gkyl_array *dxdz1 = gkyl_array_new(GKYL_DOUBLE, 9*basis.num_basis, local_ext1.volume);
   gkyl_eval_on_nodes *proj1 = gkyl_eval_on_nodes_new(&grid1, &basis, 9, &proj_one, 0);
   gkyl_eval_on_nodes_advance(proj1, 0.0, &local1, dxdz1);
   gkyl_eval_on_nodes_release(proj1);
-  gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
+  // gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
 
 
 
@@ -211,7 +211,7 @@ void test_cartesian_2x_x()
   int edge2 = 0; //lower edge
   int dir = 0; // first direction
   gkyl_bc_block_tensor_advance(bt, dir, edge1, edge2, dxdz1, dzdx2, &local1, &local2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0, bt->tensor, "tji.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0, bt->tensor, "tji.gkyl");
 
   gkyl_array_release(dxdz1);
   gkyl_array_release(dzdx2);
@@ -252,14 +252,14 @@ void test_cyl_cart_2x_z()
   gkyl_eval_on_nodes *proj2 = gkyl_eval_on_nodes_new(&grid2, &basis, 9, &proj_one, 0);
   gkyl_eval_on_nodes_advance(proj2, 0.0, &local2, dzdx2);
   gkyl_eval_on_nodes_release(proj2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0, dzdx2, "dzdx2.gkyl");
 
   // Block 1 tangents
   struct gkyl_array *dxdz1 = gkyl_array_new(GKYL_DOUBLE, 9*basis.num_basis, local_ext1.volume);
   gkyl_eval_on_nodes *proj1 = gkyl_eval_on_nodes_new(&grid1, &basis, 9, &proj_cyl, 0);
   gkyl_eval_on_nodes_advance(proj1, 0.0, &local1, dxdz1);
   gkyl_eval_on_nodes_release(proj1);
-  gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
+  // gkyl_grid_sub_array_write(&grid1, &local1, 0, dxdz1, "dxdz1.gkyl");
 
 
 
@@ -268,7 +268,7 @@ void test_cyl_cart_2x_z()
   int edge2 = 1; // upper edge
   int dir = 1; // second direction
   gkyl_bc_block_tensor_advance(bt, dir, edge1, edge2, dxdz1, dzdx2, &local1, &local2);
-  gkyl_grid_sub_array_write(&grid2, &local2, 0, bt->tensor, "tji.gkyl");
+  // gkyl_grid_sub_array_write(&grid2, &local2, 0, bt->tensor, "tji.gkyl");
 
   gkyl_array_release(dxdz1);
   gkyl_array_release(dzdx2);
@@ -300,8 +300,8 @@ void test_cartesian_3x_onecell()
   gkyl_eval_on_nodes_advance(proj, 0.0, &local, dxdz);
   gkyl_eval_on_nodes_advance(proj, 0.0, &local, dzdx);
   gkyl_eval_on_nodes_release(proj);
-  gkyl_grid_sub_array_write(&grid, &local, 0, dxdz, "dxdz.gkyl");
-  gkyl_grid_sub_array_write(&grid, &local, 0, dzdx, "dzdx.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, dxdz, "dxdz.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, dzdx, "dzdx.gkyl");
 
   struct bc_block_tensor *bt = gkyl_bc_block_tensor_new(&grid, &local, &local_ext, &basis, false);
 
diff --git a/gyrokinetic/unit/ctest_boundary_flux_gyrokinetic.c b/gyrokinetic/unit/ctest_boundary_flux_gyrokinetic.c
new file mode 100644
index 0000000000..b58560fa7a
--- /dev/null
+++ b/gyrokinetic/unit/ctest_boundary_flux_gyrokinetic.c
@@ -0,0 +1,80 @@
+// Test construction of the gyrokinetic boundary-flux updater
+// (gkyl_boundary_flux_new). Verifies that the stored direction, edge,
+// number of equations, copied grid, and skin/ghost range volumes match the
+// inputs, and that the equation object is acquired (so it survives release of
+// the caller's reference). A lightweight gyrokinetic-diffusion DG equation is
+// used as the equation object.
+#include <acutest.h>
+
+#include <gkyl_basis.h>
+#include <gkyl_boundary_flux.h>
+#include <gkyl_boundary_flux_priv.h>
+#include <gkyl_dg_eqn.h>
+#include <gkyl_dg_diffusion_gyrokinetic.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_decomp.h>
+#include <gkyl_rect_grid.h>
+
+static struct gkyl_dg_eqn*
+make_diffusion_eqn(struct gkyl_range *diff_range)
+{
+  int cdim = 1, vdim = 1, poly_order = 1;
+  struct gkyl_basis basis, cbasis;
+  gkyl_cart_modal_serendip(&cbasis, cdim, poly_order);
+  gkyl_cart_modal_gkhybrid(&basis, cdim, vdim);
+
+  gkyl_range_init(diff_range, cdim, (int[]) { 0 }, (int[]) { 7 });
+  bool dir[] = { true };
+  return gkyl_dg_diffusion_gyrokinetic_new(&basis, &cbasis, true, dir, 2, diff_range, false);
+}
+
+static void
+check_boundary_flux(int dir, enum gkyl_edge_loc edge)
+{
+  // Phase-space grid: 1x1v.
+  int cells[] = { 8, 8 };
+  int ghost[] = { 1, 0 };
+  double lower[] = { 0.0, -1.0 }, upper[] = { 1.0, 1.0 };
+
+  struct gkyl_rect_grid grid;
+  struct gkyl_range local, local_ext;
+  gkyl_rect_grid_init(&grid, 2, lower, upper, cells);
+  gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
+
+  struct gkyl_range skin_r, ghost_r;
+  gkyl_skin_ghost_ranges(&skin_r, &ghost_r, dir, edge, &local_ext, ghost);
+
+  struct gkyl_range diff_range;
+  struct gkyl_dg_eqn *eqn = make_diffusion_eqn(&diff_range);
+
+  const struct gkyl_dg_eqn *eqns[] = { eqn };
+  struct gkyl_boundary_flux *bf = gkyl_boundary_flux_new(dir, edge, &grid,
+    &skin_r, &ghost_r, 1, eqns, false);
+
+  TEST_CHECK( bf != NULL );
+  TEST_CHECK( bf->dir == dir );
+  TEST_CHECK( bf->edge == edge );
+  TEST_CHECK( bf->num_eqns == 1 );
+  TEST_CHECK( bf->use_gpu == false );
+  TEST_CHECK( bf->grid.ndim == 2 );
+  TEST_CHECK( bf->grid.cells[0] == 8 );
+  TEST_CHECK( bf->skin_r.volume == skin_r.volume );
+  TEST_CHECK( bf->ghost_r.volume == ghost_r.volume );
+  TEST_CHECK( bf->eqns[0] != NULL );
+
+  // The updater acquired its own reference; release ours and confirm the
+  // equation object is still alive (num_equations readable).
+  gkyl_dg_eqn_release(eqn);
+  TEST_CHECK( bf->eqns[0]->num_equations == 1 );
+
+  gkyl_boundary_flux_release(bf);
+}
+
+void test_boundary_flux_lower() { check_boundary_flux(0, GKYL_LOWER_EDGE); }
+void test_boundary_flux_upper() { check_boundary_flux(0, GKYL_UPPER_EDGE); }
+
+TEST_LIST = {
+  { "boundary_flux_lower", test_boundary_flux_lower },
+  { "boundary_flux_upper", test_boundary_flux_upper },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_dg_calc_gk_rad_vars_drag.c b/gyrokinetic/unit/ctest_dg_calc_gk_rad_vars_drag.c
new file mode 100644
index 0000000000..dce4b15992
--- /dev/null
+++ b/gyrokinetic/unit/ctest_dg_calc_gk_rad_vars_drag.c
@@ -0,0 +1,73 @@
+// Test the gyrokinetic radiation drag-coefficient container allocator
+// (gkyl_dg_calc_gk_rad_vars_drag_new / _release). Verifies that the returned
+// per-collision structs record the correct number of densities, that each
+// per-density drag array has the requested number of components and size, and
+// that the data arrays are writable host arrays. Release must not crash.
+#include <acutest.h>
+
+#include <gkyl_array.h>
+#include <gkyl_array_ops.h>
+#include <gkyl_dg_calc_gk_rad_vars.h>
+
+void
+test_rad_drag_alloc()
+{
+  int num_collisions = 3;
+  int num_densities[] = { 1, 2, 4 };
+  int ncomp = 6;
+  long sz = 12;
+
+  struct gkyl_gk_rad_drag *drag = gkyl_dg_calc_gk_rad_vars_drag_new(
+    num_collisions, num_densities, ncomp, sz, false);
+
+  TEST_CHECK( drag != NULL );
+
+  for (int i=0; i<num_collisions; i++) {
+    TEST_CHECK( drag[i].num_dens == num_densities[i] );
+    TEST_CHECK( drag[i].data != NULL );
+    TEST_CHECK( drag[i].on_dev != NULL );
+
+    for (int n=0; n<num_densities[i]; n++) {
+      struct gkyl_array *arr = drag[i].data[n].arr;
+      TEST_CHECK( arr != NULL );
+      TEST_CHECK( arr->ncomp == (unsigned)ncomp );
+      TEST_CHECK( arr->size == (size_t)sz );
+
+      // Array should be a usable host array: write and read back.
+      gkyl_array_clear(arr, 0.0);
+      double *d = gkyl_array_fetch(arr, 0);
+      d[0] = 3.25;
+      const double *dc = gkyl_array_cfetch(arr, 0);
+      TEST_CHECK( dc[0] == 3.25 );
+    }
+  }
+
+  gkyl_dg_calc_gk_rad_vars_drag_release(drag, num_collisions, false);
+}
+
+void
+test_rad_drag_alloc_single()
+{
+  int num_collisions = 1;
+  int num_densities[] = { 5 };
+  int ncomp = 3;
+  long sz = 8;
+
+  struct gkyl_gk_rad_drag *drag = gkyl_dg_calc_gk_rad_vars_drag_new(
+    num_collisions, num_densities, ncomp, sz, false);
+
+  TEST_CHECK( drag != NULL );
+  TEST_CHECK( drag[0].num_dens == 5 );
+  for (int n=0; n<5; n++) {
+    TEST_CHECK( drag[0].data[n].arr->ncomp == (unsigned)ncomp );
+    TEST_CHECK( drag[0].data[n].arr->size == (size_t)sz );
+  }
+
+  gkyl_dg_calc_gk_rad_vars_drag_release(drag, num_collisions, false);
+}
+
+TEST_LIST = {
+  { "rad_drag_alloc", test_rad_drag_alloc },
+  { "rad_drag_alloc_single", test_rad_drag_alloc_single },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_dg_diffusion_gyrokinetic.c b/gyrokinetic/unit/ctest_dg_diffusion_gyrokinetic.c
new file mode 100644
index 0000000000..e5bb65f9c0
--- /dev/null
+++ b/gyrokinetic/unit/ctest_dg_diffusion_gyrokinetic.c
@@ -0,0 +1,127 @@
+// Test construction of the gyrokinetic diffusion DG equation object.
+// Verifies that gkyl_dg_diffusion_gyrokinetic_new returns an object with the
+// expected base-class field (num_equations), correctly recorded const-coeff
+// flag, per-direction diffusion flags and number of phase-space basis
+// functions, and that reference counting via acquire/release works.
+#include <acutest.h>
+
+#include <gkyl_basis.h>
+#include <gkyl_dg_eqn.h>
+#include <gkyl_dg_diffusion_gyrokinetic.h>
+#include <gkyl_dg_diffusion_gyrokinetic_priv.h>
+#include <gkyl_range.h>
+
+static void
+check_diffusion(int cdim, int vdim, int poly_order, bool const_coeff, const bool *diff_in_dir)
+{
+  int pdim = cdim + vdim;
+
+  struct gkyl_basis basis, cbasis;
+  gkyl_cart_modal_serendip(&cbasis, cdim, poly_order);
+  if (poly_order == 1)
+    gkyl_cart_modal_gkhybrid(&basis, cdim, vdim);
+  else
+    gkyl_cart_modal_serendip(&basis, pdim, poly_order);
+
+  // Conf-space range used to index the diffusion coefficient.
+  struct gkyl_range diff_range;
+  int lower[GKYL_MAX_CDIM], upper[GKYL_MAX_CDIM];
+  for (int d=0; d<cdim; d++) { lower[d] = 0; upper[d] = 3; }
+  gkyl_range_init(&diff_range, cdim, lower, upper);
+
+  int diff_order = 2;
+  struct gkyl_dg_eqn *eqn = gkyl_dg_diffusion_gyrokinetic_new(&basis, &cbasis,
+    const_coeff, diff_in_dir, diff_order, &diff_range, false);
+
+  TEST_CHECK( eqn != NULL );
+  TEST_CHECK( eqn->num_equations == 1 );
+
+  // Inspect the derived object (test-only access).
+  struct dg_diffusion_gyrokinetic *diffusion =
+    container_of(eqn, struct dg_diffusion_gyrokinetic, eqn);
+
+  TEST_CHECK( diffusion->const_coeff == const_coeff );
+  TEST_CHECK( diffusion->num_basis == basis.num_basis );
+  for (int d=0; d<cdim; d++)
+    TEST_CHECK( diffusion->diff_in_dir[d] == diff_in_dir[d] );
+
+  // Volume and surface kernels must be wired up.
+  TEST_CHECK( eqn->vol_term != NULL );
+  TEST_CHECK( eqn->surf_term != NULL );
+  TEST_CHECK( eqn->boundary_surf_term != NULL );
+  TEST_CHECK( diffusion->surf[0] != NULL );
+
+  gkyl_dg_eqn_release(eqn);
+}
+
+void test_diff_ctor_1x1v_p1_const()
+{
+  bool dir[] = { true };
+  check_diffusion(1, 1, 1, true, dir);
+}
+
+void test_diff_ctor_1x2v_p1_var()
+{
+  bool dir[] = { true };
+  check_diffusion(1, 2, 1, false, dir);
+}
+
+void test_diff_ctor_2x2v_p1_xy()
+{
+  bool dir[] = { true, true };
+  check_diffusion(2, 2, 1, true, dir);
+}
+
+void test_diff_ctor_2x2v_p1_x_only()
+{
+  bool dir[] = { true, false };
+  check_diffusion(2, 2, 1, true, dir);
+}
+
+void test_diff_ctor_3x2v_p1()
+{
+  bool dir[] = { true, true, true };
+  check_diffusion(3, 2, 1, true, dir);
+}
+
+void test_diff_ctor_1x1v_p2_const()
+{
+  bool dir[] = { true };
+  check_diffusion(1, 1, 2, true, dir);
+}
+
+void
+test_diff_ctor_acquire()
+{
+  struct gkyl_basis basis, cbasis;
+  gkyl_cart_modal_serendip(&cbasis, 1, 1);
+  gkyl_cart_modal_gkhybrid(&basis, 1, 1);
+
+  struct gkyl_range diff_range;
+  gkyl_range_init(&diff_range, 1, (int[]) { 0 }, (int[]) { 3 });
+
+  bool dir[] = { true };
+  struct gkyl_dg_eqn *eqn = gkyl_dg_diffusion_gyrokinetic_new(&basis, &cbasis,
+    true, dir, 2, &diff_range, false);
+
+  struct gkyl_dg_eqn *eqn2 = gkyl_dg_eqn_acquire(eqn);
+  TEST_CHECK( eqn2 == eqn );
+  TEST_CHECK( eqn2->num_equations == 1 );
+
+  // Drop first ref; object must survive via the second reference.
+  gkyl_dg_eqn_release(eqn);
+  TEST_CHECK( eqn2->num_equations == 1 );
+
+  gkyl_dg_eqn_release(eqn2);
+}
+
+TEST_LIST = {
+  { "diff_ctor_1x1v_p1_const", test_diff_ctor_1x1v_p1_const },
+  { "diff_ctor_1x2v_p1_var", test_diff_ctor_1x2v_p1_var },
+  { "diff_ctor_2x2v_p1_xy", test_diff_ctor_2x2v_p1_xy },
+  { "diff_ctor_2x2v_p1_x_only", test_diff_ctor_2x2v_p1_x_only },
+  { "diff_ctor_3x2v_p1", test_diff_ctor_3x2v_p1 },
+  { "diff_ctor_1x1v_p2_const", test_diff_ctor_1x1v_p2_const },
+  { "diff_ctor_acquire", test_diff_ctor_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_dg_interpolate.c b/gyrokinetic/unit/ctest_dg_interpolate.c
index 03dd8a703f..0b4e984643 100644
--- a/gyrokinetic/unit/ctest_dg_interpolate.c
+++ b/gyrokinetic/unit/ctest_dg_interpolate.c
@@ -818,6 +818,32 @@ void eval_bfield_3x(double t, const double *xn, double* restrict fout, void *ctx
   fout[2] = B0;
 }
 
+// Scalar magnetic-field magnitude (|B|) evaluators. The eval_bfield_*x functions
+// above return the 3-component field vector expected by the geometry's bfield_func
+// (which takes |B| = sqrt(B.B)). When |B| is needed as a scalar field (e.g. for
+// projection with num_ret_vals=1, or inside the distribution functions), use these
+// wrappers so we don't write past a single-component output buffer.
+void eval_bmag_1x(double t, const double *xn, double* restrict fout, void *ctx)
+{
+  double B[3] = {0.0};
+  eval_bfield_1x(t, xn, B, ctx);
+  fout[0] = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
+}
+
+void eval_bmag_2x(double t, const double *xn, double* restrict fout, void *ctx)
+{
+  double B[3] = {0.0};
+  eval_bfield_2x(t, xn, B, ctx);
+  fout[0] = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
+}
+
+void eval_bmag_3x(double t, const double *xn, double* restrict fout, void *ctx)
+{
+  double B[3] = {0.0};
+  eval_bfield_3x(t, xn, B, ctx);
+  fout[0] = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
+}
+
 void eval_distf_1x1v_gk(double t, const double *xn, double* restrict fout, void *ctx)
 {
   double x = xn[0], vpar = xn[1];
@@ -974,6 +1000,15 @@ test_1x1v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   struct gkyl_range local, local_ext; // local, local-ext phase-space ranges
   gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
 
+  // Create bmag arrays.
+  struct gkyl_array *bmag = mkarr(use_gpu, confBasis.num_basis, confLocal_ext.volume);
+  struct gkyl_array *bmag_ho = use_gpu? mkarr(false, bmag->ncomp, bmag->size)
+                                      : gkyl_array_acquire(bmag);
+  gkyl_proj_on_basis *proj_bmag = gkyl_proj_on_basis_new(&confGrid, &confBasis,
+    poly_order+1, 1, eval_bmag_1x, &proj_ctx);
+  gkyl_proj_on_basis_advance(proj_bmag, 0.0, &confLocal, bmag_ho);
+  gkyl_array_copy(bmag, bmag_ho);
+
   // Create distribution function arrays.
   struct gkyl_array *distf = mkarr(use_gpu, basis.num_basis, local_ext.volume);
   struct gkyl_array *distf_ho = use_gpu? mkarr(false, distf->ncomp, distf->size)
@@ -1127,8 +1162,11 @@ test_1x1v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   gkyl_array_release(moms);
   gkyl_velocity_map_release(gvm);
   gkyl_gk_geometry_release(gk_geom);
+  gkyl_array_release(bmag);
   gkyl_array_release(distf);
+  gkyl_array_release(bmag_ho);
   gkyl_array_release(distf_ho);
+  gkyl_proj_on_basis_release(proj_bmag);
   gkyl_proj_on_basis_release(proj_distf);
 }
 
@@ -1155,11 +1193,10 @@ void eval_distf_1x2v_gk(double t, const double *xn, double* restrict fout, void
 
   double vtsq = temp/mass;
 
-  double bfield[3] = {0.0};
-  eval_bfield_1x(t, xn, bfield, ctx);
-  double bmag = sqrt(bfield[0]*bfield[0]+bfield[1]*bfield[1]+bfield[2]*bfield[2]);
+  double bmag[1] = {-1.0};
+  eval_bmag_1x(t, xn, bmag, ctx);
 
-  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag/mass)/(2.0*vtsq));
+  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag[0]/mass)/(2.0*vtsq));
 }
 
 void
@@ -1236,6 +1273,15 @@ test_1x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   struct gkyl_range local, local_ext; // local, local-ext phase-space ranges
   gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
 
+  // Create bmag arrays.
+  struct gkyl_array *bmag = mkarr(use_gpu, confBasis.num_basis, confLocal_ext.volume);
+  struct gkyl_array *bmag_ho = use_gpu? mkarr(false, bmag->ncomp, bmag->size)
+                                      : gkyl_array_acquire(bmag);
+  gkyl_proj_on_basis *proj_bmag = gkyl_proj_on_basis_new(&confGrid, &confBasis,
+    poly_order+1, 1, eval_bmag_1x, &proj_ctx);
+  gkyl_proj_on_basis_advance(proj_bmag, 0.0, &confLocal, bmag_ho);
+  gkyl_array_copy(bmag, bmag_ho);
+
   // Create distribution function arrays.
   struct gkyl_array *distf = mkarr(use_gpu, basis.num_basis, local_ext.volume);
   struct gkyl_array *distf_ho = use_gpu? mkarr(false, distf->ncomp, distf->size)
@@ -1386,8 +1432,11 @@ test_1x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   gkyl_array_release(moms);
   gkyl_velocity_map_release(gvm);
   gkyl_gk_geometry_release(gk_geom);
+  gkyl_array_release(bmag);
   gkyl_array_release(distf);
+  gkyl_array_release(bmag_ho);
   gkyl_array_release(distf_ho);
+  gkyl_proj_on_basis_release(proj_bmag);
   gkyl_proj_on_basis_release(proj_distf);
 }
 
@@ -1414,11 +1463,10 @@ void eval_distf_2x2v_gk(double t, const double *xn, double* restrict fout, void
 
   double vtsq = temp/mass;
 
-  double bfield[3] = {0.0};
-  eval_bfield_2x(t, xn, bfield, ctx);
-  double bmag = sqrt(bfield[0]*bfield[0]+bfield[1]*bfield[1]+bfield[2]*bfield[2]);
+  double bmag[1] = {-1.0};
+  eval_bmag_2x(t, xn, bmag, ctx);
 
-  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag/mass)/(2.0*vtsq));
+  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag[0]/mass)/(2.0*vtsq));
 }
 
 void
@@ -1504,6 +1552,15 @@ test_2x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   struct gkyl_range local, local_ext; // local, local-ext phase-space ranges
   gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
 
+  // Create bmag arrays.
+  struct gkyl_array *bmag = mkarr(use_gpu, confBasis.num_basis, confLocal_ext.volume);
+  struct gkyl_array *bmag_ho = use_gpu? mkarr(false, bmag->ncomp, bmag->size)
+                                      : gkyl_array_acquire(bmag);
+  gkyl_proj_on_basis *proj_bmag = gkyl_proj_on_basis_new(&confGrid, &confBasis,
+    poly_order+1, 1, eval_bmag_2x, &proj_ctx);
+  gkyl_proj_on_basis_advance(proj_bmag, 0.0, &confLocal, bmag_ho);
+  gkyl_array_copy(bmag, bmag_ho);
+
   // Create distribution function arrays.
   struct gkyl_array *distf = mkarr(use_gpu, basis.num_basis, local_ext.volume);
   struct gkyl_array *distf_ho = use_gpu? mkarr(false, distf->ncomp, distf->size)
@@ -1654,8 +1711,11 @@ test_2x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   gkyl_array_release(moms);
   gkyl_velocity_map_release(gvm);
   gkyl_gk_geometry_release(gk_geom);
+  gkyl_array_release(bmag);
   gkyl_array_release(distf);
+  gkyl_array_release(bmag_ho);
   gkyl_array_release(distf_ho);
+  gkyl_proj_on_basis_release(proj_bmag);
   gkyl_proj_on_basis_release(proj_distf);
 }
 
@@ -1682,11 +1742,10 @@ void eval_distf_3x2v_gk(double t, const double *xn, double* restrict fout, void
 
   double vtsq = temp/mass;
 
-  double bfield[3] = {0.0};
-  eval_bfield_3x(t, xn, bfield, ctx);
-  double bmag = sqrt(bfield[0]*bfield[0]+bfield[1]*bfield[1]+bfield[2]*bfield[2]);
+  double bmag[1] = {-1.0};
+  eval_bmag_3x(t, xn, bmag, ctx);
 
-  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag/mass)/(2.0*vtsq));
+  fout[0] = (den/pow(2.0*M_PI*vtsq,vdim/2.0)) * exp(-(pow(vpar-upar,2)+2.0*mu*bmag[0]/mass)/(2.0*vtsq));
 }
 
 void
@@ -1786,6 +1845,15 @@ test_3x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   struct gkyl_range local, local_ext; // local, local-ext phase-space ranges
   gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
 
+  // Create bmag arrays.
+  struct gkyl_array *bmag = mkarr(use_gpu, confBasis.num_basis, confLocal_ext.volume);
+  struct gkyl_array *bmag_ho = use_gpu? mkarr(false, bmag->ncomp, bmag->size)
+                                      : gkyl_array_acquire(bmag);
+  gkyl_proj_on_basis *proj_bmag = gkyl_proj_on_basis_new(&confGrid, &confBasis,
+    poly_order+1, 1, eval_bmag_3x, &proj_ctx);
+  gkyl_proj_on_basis_advance(proj_bmag, 0.0, &confLocal, bmag_ho);
+  gkyl_array_copy(bmag, bmag_ho);
+
   // Create distribution function arrays.
   struct gkyl_array *distf = mkarr(use_gpu, basis.num_basis, local_ext.volume);
   struct gkyl_array *distf_ho = use_gpu? mkarr(false, distf->ncomp, distf->size)
@@ -1936,8 +2004,11 @@ test_3x2v_gk(const int *cells, const int *cells_tar, int poly_order, bool use_gp
   gkyl_array_release(moms);
   gkyl_velocity_map_release(gvm);
   gkyl_gk_geometry_release(gk_geom);
+  gkyl_array_release(bmag);
   gkyl_array_release(distf);
+  gkyl_array_release(bmag_ho);
   gkyl_array_release(distf_ho);
+  gkyl_proj_on_basis_release(proj_bmag);
   gkyl_proj_on_basis_release(proj_distf);
 }
 
diff --git a/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_diff.c b/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_diff.c
new file mode 100644
index 0000000000..799cd217b4
--- /dev/null
+++ b/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_diff.c
@@ -0,0 +1,170 @@
+// Test construction of the gyrokinetic LBO diff DG equation object
+// (gkyl_dg_lbo_gyrokinetic_diff_new). Verifies the base num_equations, the
+// recorded cdim/pdim, species mass, number of conf basis functions, conf-range
+// volume and stored geometry/velocity-map pointers, plus that the surface and
+// volume kernels are wired up. Reuses the mapc2p geometry + identity
+// velocity-map boilerplate from ctest_dg_gyrokinetic.c.
+#include <acutest.h>
+
+#include <gkyl_alloc.h>
+#include <gkyl_basis.h>
+#include <gkyl_dg_eqn.h>
+#include <gkyl_dg_lbo_gyrokinetic_diff.h>
+#include <gkyl_dg_lbo_gyrokinetic_diff_priv.h>
+#include <gkyl_gk_geometry.h>
+#include <gkyl_gk_geometry_mapc2p.h>
+#include <gkyl_position_map.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_range.h>
+#include <gkyl_velocity_map.h>
+
+static void
+mapc2p(double t, const double *xc, double* GKYL_RESTRICT xp, void *ctx)
+{
+  xp[0] = xc[0]; xp[1] = xc[1]; xp[2] = xc[2];
+}
+
+static void
+bfield_func(double t, const double *xc, double* GKYL_RESTRICT fout, void *ctx)
+{
+  fout[0] = 0.0; fout[1] = 0.0; fout[2] = 1.0;
+}
+
+struct lbo_env {
+  struct gkyl_basis basis, confBasis;
+  struct gkyl_range confRange, confRange_ext;
+  struct gkyl_range phaseRange, phaseRange_ext;
+  struct gkyl_rect_grid confGrid, phaseGrid;
+  struct gkyl_position_map *pmap;
+  struct gk_geometry *gk_geom;
+  struct gkyl_velocity_map *gvm;
+};
+
+static void
+lbo_env_init(struct lbo_env *e)
+{
+  int cdim = 3, vdim = 2;
+  int pdim = cdim + vdim;
+
+  int cells[] = {8, 8, 8, 8, 8};
+  int ghost[] = {1, 1, 1, 0, 0};
+  double lower[] = {0., 0., 0., -1., 0.};
+  double upper[] = {1., 1., 1., 1., 1.};
+
+  gkyl_rect_grid_init(&e->confGrid, cdim, lower, upper, cells);
+  gkyl_create_grid_ranges(&e->confGrid, ghost, &e->confRange_ext, &e->confRange);
+
+  gkyl_rect_grid_init(&e->phaseGrid, pdim, lower, upper, cells);
+  gkyl_create_grid_ranges(&e->phaseGrid, ghost, &e->phaseRange_ext, &e->phaseRange);
+
+  double velLower[vdim], velUpper[vdim];
+  int velCells[vdim];
+  for (int d=0; d<vdim; d++) {
+    velLower[d] = lower[cdim+d];
+    velUpper[d] = upper[cdim+d];
+    velCells[d] = cells[cdim+d];
+  }
+  struct gkyl_rect_grid velGrid;
+  int velGhost[] = { 0, 0 };
+  struct gkyl_range velLocal, velLocal_ext;
+  gkyl_rect_grid_init(&velGrid, vdim, velLower, velUpper, velCells);
+  gkyl_create_grid_ranges(&velGrid, velGhost, &velLocal_ext, &velLocal);
+
+  int poly_order = 1;
+  gkyl_cart_modal_gkhybrid(&e->basis, cdim, vdim);
+  gkyl_cart_modal_serendip(&e->confBasis, cdim, poly_order);
+
+  e->pmap = gkyl_position_map_null_new();
+
+  struct gkyl_gk_geometry_inp geometry_input = {
+    .geometry_id = GKYL_GEOMETRY_MAPC2P,
+    .world = {0.0, 0.0},
+    .mapc2p = mapc2p,
+    .c2p_ctx = 0,
+    .bfield_func = bfield_func,
+    .bfield_ctx = 0,
+    .position_map = e->pmap,
+    .grid = e->confGrid,
+    .local = e->confRange,
+    .local_ext = e->confRange_ext,
+    .global = e->confRange,
+    .global_ext = e->confRange_ext,
+    .basis = e->confBasis,
+    .geo_grid = e->confGrid,
+    .geo_local = e->confRange,
+    .geo_local_ext = e->confRange_ext,
+    .geo_global = e->confRange,
+    .geo_global_ext = e->confRange_ext,
+    .geo_basis = e->confBasis,
+  };
+  e->gk_geom = gkyl_gk_geometry_mapc2p_new(&geometry_input);
+
+  struct gkyl_mapc2p_inp c2p_in = { };
+  e->gvm = gkyl_velocity_map_new(c2p_in, e->phaseGrid, velGrid,
+    e->phaseRange, e->phaseRange_ext, velLocal, velLocal_ext, false);
+}
+
+static void
+lbo_env_release(struct lbo_env *e)
+{
+  gkyl_gk_geometry_release(e->gk_geom);
+  gkyl_position_map_release(e->pmap);
+  gkyl_velocity_map_release(e->gvm);
+}
+
+void
+test_lbo_diff_ctor()
+{
+  struct lbo_env e;
+  lbo_env_init(&e);
+
+  double mass = 2.5;
+  struct gkyl_dg_eqn *eqn = gkyl_dg_lbo_gyrokinetic_diff_new(&e.confBasis, &e.basis,
+    &e.confRange, &e.phaseGrid, mass, e.gk_geom, e.gvm, false);
+
+  TEST_CHECK( eqn != NULL );
+  TEST_CHECK( eqn->num_equations == 1 );
+  TEST_CHECK( eqn->vol_term != NULL );
+  TEST_CHECK( eqn->surf_term != NULL );
+  TEST_CHECK( eqn->boundary_surf_term != NULL );
+
+  struct dg_lbo_gyrokinetic_diff *diff =
+    container_of(eqn, struct dg_lbo_gyrokinetic_diff, eqn);
+
+  TEST_CHECK( diff->cdim == 3 );
+  TEST_CHECK( diff->pdim == 5 );
+  TEST_CHECK( diff->mass == mass );
+  TEST_CHECK( diff->num_cbasis == e.confBasis.num_basis );
+  TEST_CHECK( diff->conf_range.volume == 512 );
+  TEST_CHECK( diff->gk_geom != NULL );
+  TEST_CHECK( diff->vel_map != NULL );
+  TEST_CHECK( diff->surf[0] != NULL );
+  TEST_CHECK( diff->boundary_surf[0] != NULL );
+
+  gkyl_dg_eqn_release(eqn);
+  lbo_env_release(&e);
+}
+
+void
+test_lbo_diff_acquire()
+{
+  struct lbo_env e;
+  lbo_env_init(&e);
+
+  struct gkyl_dg_eqn *eqn = gkyl_dg_lbo_gyrokinetic_diff_new(&e.confBasis, &e.basis,
+    &e.confRange, &e.phaseGrid, 1.0, e.gk_geom, e.gvm, false);
+
+  struct gkyl_dg_eqn *eqn2 = gkyl_dg_eqn_acquire(eqn);
+  TEST_CHECK( eqn2 == eqn );
+  gkyl_dg_eqn_release(eqn);
+  TEST_CHECK( eqn2->num_equations == 1 );
+  gkyl_dg_eqn_release(eqn2);
+
+  lbo_env_release(&e);
+}
+
+TEST_LIST = {
+  { "lbo_diff_ctor", test_lbo_diff_ctor },
+  { "lbo_diff_acquire", test_lbo_diff_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_drag.c b/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_drag.c
new file mode 100644
index 0000000000..8f4a4df315
--- /dev/null
+++ b/gyrokinetic/unit/ctest_dg_lbo_gyrokinetic_drag.c
@@ -0,0 +1,170 @@
+// Test construction of the gyrokinetic LBO drag DG equation object
+// (gkyl_dg_lbo_gyrokinetic_drag_new). Verifies the base num_equations, the
+// recorded cdim/pdim, species mass, number of conf basis functions, conf-range
+// volume and stored geometry/velocity-map pointers, plus that the surface and
+// volume kernels are wired up. Reuses the mapc2p geometry + identity
+// velocity-map boilerplate from ctest_dg_gyrokinetic.c.
+#include <acutest.h>
+
+#include <gkyl_alloc.h>
+#include <gkyl_basis.h>
+#include <gkyl_dg_eqn.h>
+#include <gkyl_dg_lbo_gyrokinetic_drag.h>
+#include <gkyl_dg_lbo_gyrokinetic_drag_priv.h>
+#include <gkyl_gk_geometry.h>
+#include <gkyl_gk_geometry_mapc2p.h>
+#include <gkyl_position_map.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_range.h>
+#include <gkyl_velocity_map.h>
+
+static void
+mapc2p(double t, const double *xc, double* GKYL_RESTRICT xp, void *ctx)
+{
+  xp[0] = xc[0]; xp[1] = xc[1]; xp[2] = xc[2];
+}
+
+static void
+bfield_func(double t, const double *xc, double* GKYL_RESTRICT fout, void *ctx)
+{
+  fout[0] = 0.0; fout[1] = 0.0; fout[2] = 1.0;
+}
+
+struct lbo_env {
+  struct gkyl_basis basis, confBasis;
+  struct gkyl_range confRange, confRange_ext;
+  struct gkyl_range phaseRange, phaseRange_ext;
+  struct gkyl_rect_grid confGrid, phaseGrid;
+  struct gkyl_position_map *pmap;
+  struct gk_geometry *gk_geom;
+  struct gkyl_velocity_map *gvm;
+};
+
+static void
+lbo_env_init(struct lbo_env *e)
+{
+  int cdim = 3, vdim = 2;
+  int pdim = cdim + vdim;
+
+  int cells[] = {8, 8, 8, 8, 8};
+  int ghost[] = {1, 1, 1, 0, 0};
+  double lower[] = {0., 0., 0., -1., 0.};
+  double upper[] = {1., 1., 1., 1., 1.};
+
+  gkyl_rect_grid_init(&e->confGrid, cdim, lower, upper, cells);
+  gkyl_create_grid_ranges(&e->confGrid, ghost, &e->confRange_ext, &e->confRange);
+
+  gkyl_rect_grid_init(&e->phaseGrid, pdim, lower, upper, cells);
+  gkyl_create_grid_ranges(&e->phaseGrid, ghost, &e->phaseRange_ext, &e->phaseRange);
+
+  double velLower[vdim], velUpper[vdim];
+  int velCells[vdim];
+  for (int d=0; d<vdim; d++) {
+    velLower[d] = lower[cdim+d];
+    velUpper[d] = upper[cdim+d];
+    velCells[d] = cells[cdim+d];
+  }
+  struct gkyl_rect_grid velGrid;
+  int velGhost[] = { 0, 0 };
+  struct gkyl_range velLocal, velLocal_ext;
+  gkyl_rect_grid_init(&velGrid, vdim, velLower, velUpper, velCells);
+  gkyl_create_grid_ranges(&velGrid, velGhost, &velLocal_ext, &velLocal);
+
+  int poly_order = 1;
+  gkyl_cart_modal_gkhybrid(&e->basis, cdim, vdim);
+  gkyl_cart_modal_serendip(&e->confBasis, cdim, poly_order);
+
+  e->pmap = gkyl_position_map_null_new();
+
+  struct gkyl_gk_geometry_inp geometry_input = {
+    .geometry_id = GKYL_GEOMETRY_MAPC2P,
+    .world = {0.0, 0.0},
+    .mapc2p = mapc2p,
+    .c2p_ctx = 0,
+    .bfield_func = bfield_func,
+    .bfield_ctx = 0,
+    .position_map = e->pmap,
+    .grid = e->confGrid,
+    .local = e->confRange,
+    .local_ext = e->confRange_ext,
+    .global = e->confRange,
+    .global_ext = e->confRange_ext,
+    .basis = e->confBasis,
+    .geo_grid = e->confGrid,
+    .geo_local = e->confRange,
+    .geo_local_ext = e->confRange_ext,
+    .geo_global = e->confRange,
+    .geo_global_ext = e->confRange_ext,
+    .geo_basis = e->confBasis,
+  };
+  e->gk_geom = gkyl_gk_geometry_mapc2p_new(&geometry_input);
+
+  struct gkyl_mapc2p_inp c2p_in = { };
+  e->gvm = gkyl_velocity_map_new(c2p_in, e->phaseGrid, velGrid,
+    e->phaseRange, e->phaseRange_ext, velLocal, velLocal_ext, false);
+}
+
+static void
+lbo_env_release(struct lbo_env *e)
+{
+  gkyl_gk_geometry_release(e->gk_geom);
+  gkyl_position_map_release(e->pmap);
+  gkyl_velocity_map_release(e->gvm);
+}
+
+void
+test_lbo_drag_ctor()
+{
+  struct lbo_env e;
+  lbo_env_init(&e);
+
+  double mass = 2.5;
+  struct gkyl_dg_eqn *eqn = gkyl_dg_lbo_gyrokinetic_drag_new(&e.confBasis, &e.basis,
+    &e.confRange, &e.phaseGrid, mass, e.gk_geom, e.gvm, false);
+
+  TEST_CHECK( eqn != NULL );
+  TEST_CHECK( eqn->num_equations == 1 );
+  TEST_CHECK( eqn->vol_term != NULL );
+  TEST_CHECK( eqn->surf_term != NULL );
+  TEST_CHECK( eqn->boundary_surf_term != NULL );
+
+  struct dg_lbo_gyrokinetic_drag *drag =
+    container_of(eqn, struct dg_lbo_gyrokinetic_drag, eqn);
+
+  TEST_CHECK( drag->cdim == 3 );
+  TEST_CHECK( drag->pdim == 5 );
+  TEST_CHECK( drag->mass == mass );
+  TEST_CHECK( drag->num_cbasis == e.confBasis.num_basis );
+  TEST_CHECK( drag->conf_range.volume == 512 );
+  TEST_CHECK( drag->gk_geom != NULL );
+  TEST_CHECK( drag->vel_map != NULL );
+  TEST_CHECK( drag->surf[0] != NULL );
+  TEST_CHECK( drag->boundary_surf[0] != NULL );
+
+  gkyl_dg_eqn_release(eqn);
+  lbo_env_release(&e);
+}
+
+void
+test_lbo_drag_acquire()
+{
+  struct lbo_env e;
+  lbo_env_init(&e);
+
+  struct gkyl_dg_eqn *eqn = gkyl_dg_lbo_gyrokinetic_drag_new(&e.confBasis, &e.basis,
+    &e.confRange, &e.phaseGrid, 1.0, e.gk_geom, e.gvm, false);
+
+  struct gkyl_dg_eqn *eqn2 = gkyl_dg_eqn_acquire(eqn);
+  TEST_CHECK( eqn2 == eqn );
+  gkyl_dg_eqn_release(eqn);
+  TEST_CHECK( eqn2->num_equations == 1 );
+  gkyl_dg_eqn_release(eqn2);
+
+  lbo_env_release(&e);
+}
+
+TEST_LIST = {
+  { "lbo_drag_ctor", test_lbo_drag_ctor },
+  { "lbo_drag_acquire", test_lbo_drag_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_dg_rad_gyrokinetic.c b/gyrokinetic/unit/ctest_dg_rad_gyrokinetic.c
index 0bc1a7fc84..e9ec474d4b 100644
--- a/gyrokinetic/unit/ctest_dg_rad_gyrokinetic.c
+++ b/gyrokinetic/unit/ctest_dg_rad_gyrokinetic.c
@@ -709,10 +709,10 @@ test_2x(int poly_order, bool use_gpu, double te)
 
   gkyl_dg_updater_rad_gyrokinetic_advance(slvr, &local, f, cflrate, rhs);
   
-  gkyl_grid_sub_array_write(&grid, &local, 0, rhs, "ctest_dg_rad_gyrokinetic_2x_rhs.gkyl");
-  gkyl_grid_sub_array_write(&grid, &local, 0, nvnu, "ctest_dg_rad_gyrokinetic_2x_nvnu.gkyl");
-  gkyl_grid_sub_array_write(&grid, &local, 0, nvsqnu, "ctest_dg_rad_gyrokinetic_2x_nvsqnu.gkyl");
-  gkyl_grid_sub_array_write(&grid, &local, 0, f, "ctest_dg_rad_gyrokinetic_2x_f.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, rhs, "ctest_dg_rad_gyrokinetic_2x_rhs.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, nvnu, "ctest_dg_rad_gyrokinetic_2x_nvnu.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, nvsqnu, "ctest_dg_rad_gyrokinetic_2x_nvsqnu.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, f, "ctest_dg_rad_gyrokinetic_2x_f.gkyl");
   // Take 2nd moment of rhs to find energy loss on host
   struct gkyl_dg_updater_moment *m2_calc = gkyl_dg_updater_moment_gyrokinetic_new(&grid, &confBasis, &basis,
     &confLocal, GKYL_ELECTRON_MASS, -GKYL_ELEMENTARY_CHARGE, gvm, gk_geom, NULL, GKYL_F_MOMENT_M2, false, use_gpu);
diff --git a/gyrokinetic/unit/ctest_efit.c b/gyrokinetic/unit/ctest_efit.c
index fcfa2750d2..5ad640a986 100644
--- a/gyrokinetic/unit/ctest_efit.c
+++ b/gyrokinetic/unit/ctest_efit.c
@@ -15,6 +15,7 @@
 #include <gkyl_efit.h>
 
 void test_solovev(){
+  bool write_files = false;
 
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/solovev.geqdsk",
@@ -24,16 +25,20 @@ void test_solovev(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "solovev_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "solovev_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "solovev_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "solovev_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "solovev_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "solovev_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_step(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/step.geqdsk",
     .rz_poly_order = 2,
@@ -42,16 +47,20 @@ void test_step(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "step_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "step_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "step_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "step_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "step_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "step_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_nstxu(){
+  bool write_files = false;
+
   // Uses DN configuration by default, but one can switch to SN by changing the filepath if desired.
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/nstxu_DN.geqdsk",
@@ -61,16 +70,20 @@ void test_nstxu(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  // printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g psisep=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry, efit->psisep);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "nstxu_DN_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "nstxu_DN_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "nstxu_DN_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g psisep=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry, efit->psisep);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "nstxu_DN_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "nstxu_DN_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "nstxu_DN_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_asdex(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/asdex.geqdsk",
     .rz_poly_order = 2,
@@ -78,17 +91,20 @@ void test_asdex(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "asdex_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "asdex_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "asdex_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "asdex_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "asdex_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "asdex_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_cerfon(){
-  
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/cerfon.geqdsk",
     .rz_poly_order = 2,
@@ -97,16 +113,20 @@ void test_cerfon(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "cerfon_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "cerfon_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "cerfon_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "cerfon_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "cerfon_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "cerfon_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_elliptical(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/elliptical.geqdsk",
     .rz_poly_order = 2,
@@ -115,16 +135,20 @@ void test_elliptical(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "elliptical_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "elliptical_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "elliptical_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "elliptical_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "elliptical_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "elliptical_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_wham(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/wham.geqdsk",
     .rz_poly_order = 2,
@@ -133,10 +157,12 @@ void test_wham(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "wham_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "wham_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "wham_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "wham_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "wham_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "wham_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
@@ -144,6 +170,8 @@ void test_wham(){
 
 
 void test_tcv(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/tcv.geqdsk",
     .rz_poly_order = 2,
@@ -151,16 +179,20 @@ void test_tcv(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "tcv_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "tcv_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "tcv_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "tcv_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "tcv_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "tcv_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_mast(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/mast.geqdsk",
     .rz_poly_order = 2,
@@ -169,16 +201,20 @@ void test_mast(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "mast_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "mast_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "mast_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "mast_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "mast_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "mast_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
 }
 
 void test_ltx(){
+  bool write_files = false;
+
   struct gkyl_efit_inp inp  = {
     .filepath = "gyrokinetic/data/eqdsk/LTX_103955_03.eqdsk",
     .rz_poly_order = 2,
@@ -187,10 +223,12 @@ void test_ltx(){
   };
   struct gkyl_efit* efit = gkyl_efit_new(&inp);
 
-  //printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "ltx_psi.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "ltx_fpol.gkyl");
-  gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "ltx_q.gkyl");
+  if (write_files) {
+    printf( "rdim=%g zdim=%g rcentr=%g rleft=%g zmid=%g  rmaxis=%g zmaxis=%g simag=%1.16e sibry=%1.16e bcentr=%g  current=%g simag=%g rmaxis=%g   zmaxis=%g sibry=%g \n", efit->rdim, efit->zdim, efit->rcentr, efit->rleft, efit->zmid, efit->rmaxis, efit->zmaxis, efit->simag, efit->sibry, efit->bcentr, efit-> current, efit->simag, efit->rmaxis, efit-> zmaxis, efit->sibry);
+    gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "ltx_psi.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->fpolflux, "ltx_fpol.gkyl");
+    gkyl_grid_sub_array_write(&efit->fluxgrid, &efit->fluxlocal, 0, efit->qflux, "ltx_q.gkyl");
+  }
 
   gkyl_efit_release(efit);
 
diff --git a/gyrokinetic/unit/ctest_gk_anomalous_diffusion.c b/gyrokinetic/unit/ctest_gk_anomalous_diffusion.c
new file mode 100644
index 0000000000..bfb8991e4d
--- /dev/null
+++ b/gyrokinetic/unit/ctest_gk_anomalous_diffusion.c
@@ -0,0 +1,113 @@
+// Test construction of the gyrokinetic anomalous-diffusion DG equation object
+// (gkyl_gk_anomalous_diffusion_new). Verifies the base num_equations, the
+// number of phase-space basis functions recorded, the stored conf-space range
+// volume, and that the volume / surface / boundary kernels are wired up for
+// the supported 2x2v p1 case under a couple of boundary-condition choices.
+// Also checks acquire/release reference counting.
+#include <acutest.h>
+
+#include <gkyl_basis.h>
+#include <gkyl_dg_eqn.h>
+#include <gkyl_gk_anomalous_diffusion.h>
+#include <gkyl_gk_anomalous_diffusion_priv.h>
+#include <gkyl_gk_bc_type.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_decomp.h>
+#include <gkyl_rect_grid.h>
+
+static void
+make_ranges(int cdim, struct gkyl_range *confRange, struct gkyl_range *confRange_ext,
+  struct gkyl_basis *basis, struct gkyl_basis *cbasis)
+{
+  int vdim = 2;
+  int pdim = cdim + vdim;
+  int cells[GKYL_MAX_DIM];
+  int ghost[GKYL_MAX_DIM];
+  double lower[GKYL_MAX_DIM], upper[GKYL_MAX_DIM];
+  for (int d=0; d<cdim; d++) { cells[d] = 8; ghost[d] = 1; lower[d] = 0.0; upper[d] = 1.0; }
+
+  struct gkyl_rect_grid confGrid;
+  gkyl_rect_grid_init(&confGrid, cdim, lower, upper, cells);
+  gkyl_create_grid_ranges(&confGrid, ghost, confRange_ext, confRange);
+
+  gkyl_cart_modal_serendip(cbasis, cdim, 1);
+  gkyl_cart_modal_gkhybrid(basis, cdim, vdim);
+}
+
+void
+test_anom_diff_ctor_local()
+{
+  int cdim = 2;
+  struct gkyl_basis basis, cbasis;
+  struct gkyl_range confRange, confRange_ext;
+  make_ranges(cdim, &confRange, &confRange_ext, &basis, &cbasis);
+
+  // "ELSE" branch (local boundary kernels) using copy BCs.
+  struct gkyl_dg_eqn *eqn = gkyl_gk_anomalous_diffusion_new(&basis, &cbasis,
+    &confRange, GKYL_BC_GK_SPECIES_COPY, GKYL_BC_GK_SPECIES_COPY, false);
+
+  TEST_CHECK( eqn != NULL );
+  TEST_CHECK( eqn->num_equations == 1 );
+  TEST_CHECK( eqn->vol_term != NULL );
+  TEST_CHECK( eqn->surf_term != NULL );
+
+  struct gk_anomalous_diffusion *diffusion =
+    container_of(eqn, struct gk_anomalous_diffusion, eqn);
+
+  TEST_CHECK( diffusion->conf_range.volume == 64 );  // 8x8 interior cells
+  TEST_CHECK( diffusion->surf != NULL );
+  TEST_CHECK( diffusion->boundary_surf[0] != NULL );
+  TEST_CHECK( diffusion->boundary_surf[1] != NULL );
+
+  gkyl_dg_eqn_release(eqn);
+}
+
+void
+test_anom_diff_ctor_zeroflux()
+{
+  int cdim = 2;
+  struct gkyl_basis basis, cbasis;
+  struct gkyl_range confRange, confRange_ext;
+  make_ranges(cdim, &confRange, &confRange_ext, &basis, &cbasis);
+
+  // ZERO_FLUX branch on both ends.
+  struct gkyl_dg_eqn *eqn = gkyl_gk_anomalous_diffusion_new(&basis, &cbasis,
+    &confRange, GKYL_BC_GK_SPECIES_ZERO_FLUX, GKYL_BC_GK_SPECIES_ZERO_FLUX, false);
+
+  TEST_CHECK( eqn != NULL );
+  TEST_CHECK( eqn->num_equations == 1 );
+  TEST_CHECK( eqn->vol_term != NULL );
+  TEST_CHECK( eqn->surf_term != NULL );
+
+  struct gk_anomalous_diffusion *diffusion =
+    container_of(eqn, struct gk_anomalous_diffusion, eqn);
+  TEST_CHECK( diffusion->conf_range.volume == 64 );
+  TEST_CHECK( diffusion->surf != NULL );
+
+  gkyl_dg_eqn_release(eqn);
+}
+
+void
+test_anom_diff_acquire()
+{
+  int cdim = 2;
+  struct gkyl_basis basis, cbasis;
+  struct gkyl_range confRange, confRange_ext;
+  make_ranges(cdim, &confRange, &confRange_ext, &basis, &cbasis);
+
+  struct gkyl_dg_eqn *eqn = gkyl_gk_anomalous_diffusion_new(&basis, &cbasis,
+    &confRange, GKYL_BC_GK_SPECIES_COPY, GKYL_BC_GK_SPECIES_COPY, false);
+
+  struct gkyl_dg_eqn *eqn2 = gkyl_dg_eqn_acquire(eqn);
+  TEST_CHECK( eqn2 == eqn );
+  gkyl_dg_eqn_release(eqn);
+  TEST_CHECK( eqn2->num_equations == 1 );
+  gkyl_dg_eqn_release(eqn2);
+}
+
+TEST_LIST = {
+  { "anom_diff_ctor_local", test_anom_diff_ctor_local },
+  { "anom_diff_ctor_zeroflux", test_anom_diff_ctor_zeroflux },
+  { "anom_diff_acquire", test_anom_diff_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_gk_geometry_mirror.c b/gyrokinetic/unit/ctest_gk_geometry_mirror.c
index c5fb213b6e..873f796af8 100644
--- a/gyrokinetic/unit/ctest_gk_geometry_mirror.c
+++ b/gyrokinetic/unit/ctest_gk_geometry_mirror.c
@@ -282,7 +282,7 @@ test_3x_p1_straight_cylinder()
 
   struct gk_geometry *gk_geom = gkyl_gk_geometry_mirror_new(&geometry_input);
 
-  write_geometry(gk_geom, grid, range, "straight_cylinder");
+  // write_geometry(gk_geom, grid, range, "straight_cylinder");
 
   // Define nodal operations
   enum { PSI_IDX, AL_IDX, TH_IDX }; // arrangement of computational coordinates
diff --git a/gyrokinetic/unit/ctest_gk_neut_fluid_prim_vars.c b/gyrokinetic/unit/ctest_gk_neut_fluid_prim_vars.c
new file mode 100644
index 0000000000..c78805a203
--- /dev/null
+++ b/gyrokinetic/unit/ctest_gk_neut_fluid_prim_vars.c
@@ -0,0 +1,241 @@
+// Test the GK neutral-fluid primitive-variable updater.
+//
+// The fluid moments are (rho, rho*ux, rho*uy, rho*uz, totalE).  For spatially
+// uniform (constant) moments the primitive variables are exact algebraic
+// combinations:
+//   udrift_i = (rho*u_i)/rho
+//   p        = (gas_gamma-1)*(E - 1/2 rho u^2)
+//   T        = p/rho * mass            [ = (gas_gamma-1)*(mass*E - 1/2 (rho u)^2)/rho ]
+//   flowE    = 1/2 rho u^2
+//   thermalE = p/(gas_gamma-1) = E - 1/2 rho u^2
+// Using is_integrated=true the updater returns the cell integral of each
+// primitive variable, i.e. (value)*(cell volume), which we check exactly.
+#include <acutest.h>
+
+#include <gkyl_array.h>
+#include <gkyl_array_ops.h>
+#include <gkyl_basis.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_rect_decomp.h>
+#include <gkyl_util.h>
+#include <gkyl_gk_neut_fluid_prim_vars.h>
+#include <math.h>
+
+// Exported by the library but (in this build) not declared in the public header.
+void gkyl_gk_neut_fluid_prim_vars_flow_energy_advance(struct gkyl_gk_neut_fluid_prim_vars *up,
+  const struct gkyl_array* moms, struct gkyl_array *out, int out_coff);
+
+// Build a constant DG field of ncomp moment-components on a 1x grid; each
+// component is filled with the constant value vals[c].
+static struct gkyl_array*
+mk_const_moms(struct gkyl_basis *basis, struct gkyl_range *range, int nmom, const double *vals)
+{
+  int nb = basis->num_basis;
+  struct gkyl_array *arr = gkyl_array_new(GKYL_DOUBLE, nmom*nb, range->volume);
+  gkyl_array_clear(arr, 0.0);
+
+  // For a constant function value v, only the 0th DG coefficient is nonzero
+  // and equals v*sqrt(2)^cdim.
+  double fac = pow(sqrt(2.0), basis->ndim);
+
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, range);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(range, iter.idx);
+    double *d = gkyl_array_fetch(arr, lidx);
+    for (int c=0; c<nmom; c++)
+      d[c*nb + 0] = vals[c]*fac;
+  }
+  return arr;
+}
+
+struct setup {
+  struct gkyl_basis basis;
+  struct gkyl_rect_grid grid;
+  struct gkyl_range local, local_ext;
+  double cell_vol;
+};
+
+static void
+make_setup(struct setup *s, int poly_order)
+{
+  double lower[] = {0.0}, upper[] = {2.0};
+  int cells[] = {4};
+  int dim = 1;
+  gkyl_rect_grid_init(&s->grid, dim, lower, upper, cells);
+  int ghost[] = {1};
+  gkyl_create_grid_ranges(&s->grid, ghost, &s->local_ext, &s->local);
+  gkyl_cart_modal_serendip(&s->basis, dim, poly_order);
+  s->cell_vol = (upper[0]-lower[0])/cells[0];
+}
+
+// Read the integrated (cell-integral) scalar from component c of out, at the
+// first cell of the local range, and divide out the cell volume to recover the
+// pointwise primitive-variable value.
+static double
+read_val(struct gkyl_array *out, struct gkyl_range *range, int c, double cell_vol)
+{
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, range);
+  gkyl_range_iter_next(&iter);
+  long lidx = gkyl_range_idx(range, iter.idx);
+  const double *d = gkyl_array_cfetch(out, lidx);
+  return d[c]/cell_vol;
+}
+
+void
+test_prim_vars_udrift_pressure_temp()
+{
+  struct setup s;
+  make_setup(&s, 1);
+
+  double gas_gamma = 5.0/3.0;
+  double mass = 2.0;
+
+  // Choose moments.
+  double rho = 3.0, ux = 1.5, uy = -0.5, uz = 0.25;
+  double rhoux = rho*ux, rhouy = rho*uy, rhouz = rho*uz;
+  double usq = ux*ux + uy*uy + uz*uz;
+  double p = 4.0; // desired pressure
+  // E = p/(gas_gamma-1) + 1/2 rho u^2.
+  double E = p/(gas_gamma-1.0) + 0.5*rho*usq;
+  double moms[] = {rho, rhoux, rhouy, rhouz, E};
+
+  struct gkyl_array *m = mk_const_moms(&s.basis, &s.local_ext, 5, moms);
+
+  double T = p/rho*mass;             // T = mass*p/rho per kernel definition
+  double flowE = 0.5*rho*usq;
+
+  // --- udrift ---
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_UDRIFT, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 3, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_udrift_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), ux, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 1, s.cell_vol), uy, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 2, s.cell_vol), uz, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  // --- pressure ---
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_PRESSURE, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 1, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_pressure_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), p, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  // --- temperature ---
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_TEMP, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 1, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_temp_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), T, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  // --- flow energy ---
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_FLOW_ENERGY, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 1, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_flow_energy_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), flowE, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  gkyl_array_release(m);
+}
+
+void
+test_prim_vars_combined()
+{
+  struct setup s;
+  make_setup(&s, 1);
+
+  double gas_gamma = 1.4;
+  double mass = 1.0;
+
+  double rho = 2.0, ux = 0.5, uy = 1.0, uz = -1.5;
+  double usq = ux*ux + uy*uy + uz*uz;
+  double p = 6.0;
+  double E = p/(gas_gamma-1.0) + 0.5*rho*usq;
+  double moms[] = {rho, rho*ux, rho*uy, rho*uz, E};
+  double T = p/rho*mass;
+
+  struct gkyl_array *m = mk_const_moms(&s.basis, &s.local_ext, 5, moms);
+
+  // udrift + pressure -> (ux, uy, uz, p)
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_UDRIFT_PRESSURE, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 4, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_udrift_pressure_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), ux, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 1, s.cell_vol), uy, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 2, s.cell_vol), uz, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 3, s.cell_vol), p, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  // udrift + temperature -> (ux, uy, uz, T)
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_UDRIFT_TEMP, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 4, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_udrift_temp_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), ux, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 1, s.cell_vol), uy, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 2, s.cell_vol), uz, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 3, s.cell_vol), T, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  // LTE -> (n=rho/mass, ux, uy, uz, T/mass)
+  {
+    struct gkyl_gk_neut_fluid_prim_vars *up = gkyl_gk_neut_fluid_prim_vars_new(
+      gas_gamma, mass, &s.basis, &s.grid, &s.local,
+      GKYL_GK_NEUT_FLUID_PRIM_VARS_LTE, true, false);
+    struct gkyl_array *out = gkyl_array_new(GKYL_DOUBLE, 5, s.local_ext.volume);
+    gkyl_array_clear(out, 0.0);
+    gkyl_gk_neut_fluid_prim_vars_lte_advance(up, m, out, 0);
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 0, s.cell_vol), rho/mass, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 1, s.cell_vol), ux, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 2, s.cell_vol), uy, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 3, s.cell_vol), uz, 1e-12) );
+    TEST_CHECK( gkyl_compare(read_val(out, &s.local, 4, s.cell_vol), T/mass, 1e-12) );
+    gkyl_array_release(out);
+    gkyl_gk_neut_fluid_prim_vars_release(up);
+  }
+
+  gkyl_array_release(m);
+}
+
+TEST_LIST = {
+  { "prim_vars_udrift_pressure_temp", test_prim_vars_udrift_pressure_temp },
+  { "prim_vars_combined", test_prim_vars_combined },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_gkgeom.c b/gyrokinetic/unit/ctest_gkgeom.c
index df401cc0c5..c1b48be4d7 100644
--- a/gyrokinetic/unit/ctest_gkgeom.c
+++ b/gyrokinetic/unit/ctest_gkgeom.c
@@ -57,7 +57,7 @@ ellip_unit(void)
   gkyl_eval_on_nodes_advance(eon, 0.0, &rzlocal, psiRZ);
   gkyl_eval_on_nodes_release(eon);
 
-  gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "ellip_psi.gkyl");
+  // gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "ellip_psi.gkyl");
 
   gkyl_gkgeom *geo = gkyl_gkgeom_new(&(struct gkyl_gkgeom_inp) {
       // psiRZ and related inputs
@@ -155,7 +155,7 @@ cerfon_unit(void)
   gkyl_eval_on_nodes_advance(eon, 0.0, &rzlocal, psiRZ);
   gkyl_eval_on_nodes_release(eon);
 
-  gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "cerfon_psi.gkyl");
+  // gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "cerfon_psi.gkyl");
 
   gkyl_gkgeom *geo = gkyl_gkgeom_new(&(struct gkyl_gkgeom_inp) {
       // psiRZ and related inputs
@@ -351,7 +351,7 @@ wham_2l_unit(void)
   gkyl_eval_on_nodes_advance(eon, 0.0, &rzlocal, psiRZ);
   gkyl_eval_on_nodes_release(eon);
 
-  gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "wham_psi.gkyl");
+  // gkyl_grid_sub_array_write(&rzgrid, &rzlocal, 0, psiRZ, "wham_psi.gkyl");
 
   gkyl_gkgeom *geo = gkyl_gkgeom_new(&(struct gkyl_gkgeom_inp) {
       // psiRZ and related inputs
diff --git a/gyrokinetic/unit/ctest_gyrokinetic_cross_prim_moms_bgk.c b/gyrokinetic/unit/ctest_gyrokinetic_cross_prim_moms_bgk.c
index f626ea5b16..da5e7d2411 100644
--- a/gyrokinetic/unit/ctest_gyrokinetic_cross_prim_moms_bgk.c
+++ b/gyrokinetic/unit/ctest_gyrokinetic_cross_prim_moms_bgk.c
@@ -458,10 +458,109 @@ void test_1x2v(int poly_order, bool use_gpu)
   gkyl_proj_on_basis_release(proj_vtsq_i);
 }
 
+// Test the gyrokinetic BGK cross primitive-moment calculator in the special
+// case of two identical species.  The cross primitive moments (n, upar, vt^2)
+// of a species relaxing against an identical species must be identical to its
+// own primitive moments: the density is unchanged, and the cross drift speed
+// and thermal speed equal the (common) input values.  Uses spatially-constant
+// fields so the expected DG coefficients are exact.
+static void
+run_equal_species(int cdim, int vdim, int poly_order,
+  double n0, double upar0, double vtsq0)
+{
+  int ndim = cdim + vdim;
+
+  double lower[GKYL_MAX_DIM], upper[GKYL_MAX_DIM];
+  int cells[GKYL_MAX_DIM];
+  for (int d=0; d<cdim; d++) { lower[d] = -0.5; upper[d] = 0.5; cells[d] = 2; }
+  for (int d=0; d<vdim; d++) {
+    lower[cdim+d] = (d==0)? -6.0 : 0.0;
+    upper[cdim+d] = (d==0)?  6.0 : 6.0;
+    cells[cdim+d] = 8;
+  }
+
+  double confLower[GKYL_MAX_CDIM], confUpper[GKYL_MAX_CDIM];
+  int confCells[GKYL_MAX_CDIM];
+  for (int d=0; d<cdim; d++) {
+    confLower[d] = lower[d]; confUpper[d] = upper[d]; confCells[d] = cells[d];
+  }
+
+  struct gkyl_rect_grid grid, confGrid;
+  gkyl_rect_grid_init(&grid, ndim, lower, upper, cells);
+  gkyl_rect_grid_init(&confGrid, cdim, confLower, confUpper, confCells);
+
+  struct gkyl_basis basis, confBasis;
+  if (poly_order == 1)
+    gkyl_cart_modal_gkhybrid(&basis, cdim, vdim);
+  else
+    gkyl_cart_modal_serendip(&basis, ndim, poly_order);
+  gkyl_cart_modal_serendip(&confBasis, cdim, poly_order);
+
+  int confGhost[GKYL_MAX_CDIM];
+  for (int d=0; d<cdim; d++) confGhost[d] = 1;
+  struct gkyl_range confLocal, confLocal_ext;
+  gkyl_create_grid_ranges(&confGrid, confGhost, &confLocal_ext, &confLocal);
+
+  int nb = confBasis.num_basis;
+  double fac = pow(sqrt(2.0), cdim); // 0th coeff scale of a constant
+
+  // prim_moms layout: [den | upar | vtsq], each nb coefficients.
+  struct gkyl_array *prim_moms = mkarr(3*nb, confLocal_ext.volume, false);
+  gkyl_array_clear(prim_moms, 0.0);
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &confLocal_ext);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(&confLocal_ext, iter.idx);
+    double *d = gkyl_array_fetch(prim_moms, lidx);
+    d[0*nb + 0] = n0*fac;
+    d[1*nb + 0] = upar0*fac;
+    d[2*nb + 0] = vtsq0*fac;
+  }
+
+  double mass = 1.5;
+  double delta_sr = 1.0, betap1 = 1.0;
+
+  gkyl_gyrokinetic_cross_prim_moms_bgk *calc =
+    gkyl_gyrokinetic_cross_prim_moms_bgk_new(&basis, &confBasis, false);
+
+  struct gkyl_array *out = mkarr(3*nb, confLocal_ext.volume, false);
+  gkyl_array_clear(out, 0.0);
+
+  // Self and other are identical.
+  gkyl_gyrokinetic_cross_prim_moms_bgk_advance(calc, &confLocal, delta_sr, betap1,
+    mass, prim_moms, mass, prim_moms, out);
+
+  gkyl_gyrokinetic_cross_prim_moms_bgk_release(calc);
+
+  gkyl_range_iter_init(&iter, &confLocal);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(&confLocal, iter.idx);
+    const double *o = gkyl_array_cfetch(out, lidx);
+    // Recover pointwise values from 0th coefficients.
+    double n_out    = o[0*nb + 0]/fac;
+    double upar_out = o[1*nb + 0]/fac;
+    double vtsq_out = o[2*nb + 0]/fac;
+    TEST_CHECK( gkyl_compare(n_out, n0, 1e-10*fabs(n0)) );
+    TEST_CHECK( gkyl_compare(upar_out, upar0, 1e-9*fmax(fabs(upar0), 1.0)) );
+    TEST_CHECK( gkyl_compare(vtsq_out, vtsq0, 1e-9*fabs(vtsq0)) );
+    TEST_MSG("n=%g upar=%g vtsq=%g", n_out, upar_out, vtsq_out);
+  }
+
+  gkyl_array_release(prim_moms);
+  gkyl_array_release(out);
+}
+
+void test_equal_1x1v_p1() { run_equal_species(1, 1, 1, 2.0e19, 1.0e4, 5.0e11); }
+void test_equal_1x2v_p1() { run_equal_species(1, 2, 1, 3.0e19, -2.0e4, 8.0e11); }
+void test_equal_1x1v_zeroflow_p1() { run_equal_species(1, 1, 1, 1.0e19, 0.0, 1.0e12); }
+
 void test_1x1v_p1() {test_1x1v(1, false);}
 void test_1x2v_p1() {test_1x2v(1, false);}
 TEST_LIST = {
   {"test_1x1v_p1", test_1x1v_p1},
   {"test_1x2v_p1", test_1x2v_p1},
+  { "equal_1x1v_p1", test_equal_1x1v_p1 },
+  { "equal_1x2v_p1", test_equal_1x2v_p1 },
+  { "equal_1x1v_zeroflow_p1", test_equal_1x1v_zeroflow_p1 },
   {NULL, NULL},
 };
diff --git a/gyrokinetic/unit/ctest_loss_cone_mask_gyrokinetic.c b/gyrokinetic/unit/ctest_loss_cone_mask_gyrokinetic.c
index fa3caa80a7..98417d60ad 100644
--- a/gyrokinetic/unit/ctest_loss_cone_mask_gyrokinetic.c
+++ b/gyrokinetic/unit/ctest_loss_cone_mask_gyrokinetic.c
@@ -376,10 +376,10 @@ test_1x2v_gk(int poly_order, bool use_gpu)
     sprintf(fname, "ctest_loss_cone_mask_gyrokinetic_1x2v_p%d_dev.gkyl", poly_order);
   else
     sprintf(fname, "ctest_loss_cone_mask_gyrokinetic_1x2v_p%d_ho.gkyl", poly_order);
-  gkyl_grid_sub_array_write(&grid, &local, 0, mask_ho, fname);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, mask_ho, fname);
 
   sprintf(fname, "ctest_loss_cone_mask_gyrokinetic_1x2v_p%d_ref.gkyl", poly_order);
-  gkyl_grid_sub_array_write(&grid, &local, 0, mask_ref_ho, fname);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, mask_ref_ho, fname);
 
   if (use_gpu) {
     gkyl_cu_free(bmag_max);
diff --git a/gyrokinetic/unit/ctest_ltx_miller.c b/gyrokinetic/unit/ctest_ltx_miller.c
index 2f680d1103..f19adf3d41 100644
--- a/gyrokinetic/unit/ctest_ltx_miller.c
+++ b/gyrokinetic/unit/ctest_ltx_miller.c
@@ -141,7 +141,7 @@ test_ltx_miller()
   };
 
   struct gk_geometry* up = gkyl_gk_geometry_tok_new(&geometry_inp); 
-  write_geometry(up, cgrid, cbasis, clocal, "ltx_miller");
+  // write_geometry(up, cgrid, cbasis, clocal, "ltx_miller");
   gkyl_gk_geometry_release(up);
   gkyl_position_map_release(pmap);
 
diff --git a/gyrokinetic/unit/ctest_mom_bcorr_lbo_gyrokinetic.c b/gyrokinetic/unit/ctest_mom_bcorr_lbo_gyrokinetic.c
new file mode 100644
index 0000000000..7069b0dc80
--- /dev/null
+++ b/gyrokinetic/unit/ctest_mom_bcorr_lbo_gyrokinetic.c
@@ -0,0 +1,135 @@
+// Test construction of the gyrokinetic LBO boundary-correction moment type.
+// Verifies the returned gkyl_mom_type carries the expected dim/basis fields,
+// that num_mom == 2 (parallel-velocity + energy corrections), and that
+// acquire/release reference counting works. Uses an identity velocity map.
+#include <acutest.h>
+
+#include <gkyl_basis.h>
+#include <gkyl_mom_type.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_velocity_map.h>
+#include <gkyl_mom_bcorr_lbo_gyrokinetic.h>
+
+static struct gkyl_velocity_map*
+make_identity_vmap(int cdim, int vdim, double *plower, double *pupper, int *pcells)
+{
+  int pdim = cdim + vdim;
+
+  double vlower[3], vupper[3];
+  int vcells[3];
+  for (int d=0; d<vdim; d++) {
+    vlower[d] = plower[cdim+d];
+    vupper[d] = pupper[cdim+d];
+    vcells[d] = pcells[cdim+d];
+  }
+
+  struct gkyl_rect_grid grid, vgrid;
+  gkyl_rect_grid_init(&grid, pdim, plower, pupper, pcells);
+  gkyl_rect_grid_init(&vgrid, vdim, vlower, vupper, vcells);
+
+  int ghost[GKYL_MAX_DIM] = {0};
+  struct gkyl_range local, local_ext;
+  gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
+
+  int vghost[3] = {0};
+  struct gkyl_range vlocal, vlocal_ext;
+  gkyl_create_grid_ranges(&vgrid, vghost, &vlocal_ext, &vlocal);
+
+  struct gkyl_mapc2p_inp c2p_in = { };
+  return gkyl_velocity_map_new(c2p_in, grid, vgrid, local, local_ext,
+    vlocal, vlocal_ext, false);
+}
+
+static void
+check_bcorr(int cdim, int vdim, int poly_order, double *plower, double *pupper, int *pcells)
+{
+  int pdim = cdim + vdim;
+
+  struct gkyl_basis cbasis, pbasis;
+  gkyl_cart_modal_serendip(&cbasis, cdim, poly_order);
+  if (poly_order == 1)
+    gkyl_cart_modal_gkhybrid(&pbasis, cdim, vdim);
+  else
+    gkyl_cart_modal_serendip(&pbasis, pdim, poly_order);
+
+  struct gkyl_velocity_map *gvm = make_identity_vmap(cdim, vdim, plower, pupper, pcells);
+
+  double mass = 1.25;
+  struct gkyl_mom_type *bcorr = gkyl_mom_bcorr_lbo_gyrokinetic_new(&cbasis, &pbasis,
+    mass, gvm, false);
+
+  TEST_CHECK( bcorr != NULL );
+  TEST_CHECK( bcorr->cdim == cdim );
+  TEST_CHECK( bcorr->pdim == pdim );
+  TEST_CHECK( bcorr->poly_order == poly_order );
+  TEST_CHECK( bcorr->num_config == cbasis.num_basis );
+  TEST_CHECK( bcorr->num_phase == pbasis.num_basis );
+  // Two boundary-correction moments: parallel velocity and energy.
+  TEST_CHECK( bcorr->num_mom == 2 );
+  TEST_CHECK( gkyl_mom_type_num_mom(bcorr) == 2 );
+
+  gkyl_mom_type_release(bcorr);
+  gkyl_velocity_map_release(gvm);
+}
+
+void test_bcorr_ctor_1x1v_p1()
+{
+  double lo[] = {-1.0, -6.0}, up[] = {1.0, 6.0};
+  int cells[] = {4, 8};
+  check_bcorr(1, 1, 1, lo, up, cells);
+}
+void test_bcorr_ctor_1x2v_p1()
+{
+  double lo[] = {-1.0, -6.0, 0.0}, up[] = {1.0, 6.0, 36.0};
+  int cells[] = {4, 8, 4};
+  check_bcorr(1, 2, 1, lo, up, cells);
+}
+void test_bcorr_ctor_2x2v_p1()
+{
+  double lo[] = {-1.0, -1.0, -6.0, 0.0}, up[] = {1.0, 1.0, 6.0, 36.0};
+  int cells[] = {4, 4, 8, 4};
+  check_bcorr(2, 2, 1, lo, up, cells);
+}
+void test_bcorr_ctor_3x2v_p1()
+{
+  double lo[] = {-1.0, -1.0, -1.0, -6.0, 0.0}, up[] = {1.0, 1.0, 1.0, 6.0, 36.0};
+  int cells[] = {2, 2, 2, 8, 4};
+  check_bcorr(3, 2, 1, lo, up, cells);
+}
+
+void
+test_bcorr_ctor_acquire()
+{
+  double lo[] = {-1.0, -6.0, 0.0}, up[] = {1.0, 6.0, 36.0};
+  int cells[] = {4, 8, 4};
+
+  struct gkyl_basis cbasis, pbasis;
+  gkyl_cart_modal_serendip(&cbasis, 1, 1);
+  gkyl_cart_modal_gkhybrid(&pbasis, 1, 2);
+
+  struct gkyl_velocity_map *gvm = make_identity_vmap(1, 2, lo, up, cells);
+
+  struct gkyl_mom_type *bcorr = gkyl_mom_bcorr_lbo_gyrokinetic_new(&cbasis, &pbasis,
+    1.0, gvm, false);
+
+  struct gkyl_mom_type *bcorr2 = gkyl_mom_type_acquire(bcorr);
+  TEST_CHECK( bcorr2 == bcorr );
+
+  gkyl_mom_type_release(bcorr);
+
+  TEST_CHECK( bcorr2->num_mom == 2 );
+  TEST_CHECK( bcorr2->cdim == 1 );
+
+  gkyl_mom_type_release(bcorr2);
+  gkyl_velocity_map_release(gvm);
+}
+
+TEST_LIST = {
+  { "bcorr_ctor_1x1v_p1", test_bcorr_ctor_1x1v_p1 },
+  { "bcorr_ctor_1x2v_p1", test_bcorr_ctor_1x2v_p1 },
+  { "bcorr_ctor_2x2v_p1", test_bcorr_ctor_2x2v_p1 },
+  { "bcorr_ctor_3x2v_p1", test_bcorr_ctor_3x2v_p1 },
+  { "bcorr_ctor_acquire", test_bcorr_ctor_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_mom_gyrokinetic_types.c b/gyrokinetic/unit/ctest_mom_gyrokinetic_types.c
new file mode 100644
index 0000000000..f9d8e19eaa
--- /dev/null
+++ b/gyrokinetic/unit/ctest_mom_gyrokinetic_types.c
@@ -0,0 +1,240 @@
+// Test construction of gyrokinetic moment-type objects for every supported
+// moment name, verifying the reported num_mom and basic dimension fields.
+// Covers both the standard (gkyl_mom_gyrokinetic_new) and integrated
+// (gkyl_int_mom_gyrokinetic_new) constructors. Host only.
+#include <acutest.h>
+
+#include <gkyl_array.h>
+#include <gkyl_basis.h>
+#include <gkyl_mom_type.h>
+#include <gkyl_mom_gyrokinetic.h>
+#include <gkyl_gk_geometry.h>
+#include <gkyl_gk_geometry_mapc2p.h>
+#include <gkyl_position_map.h>
+#include <gkyl_velocity_map.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_grid.h>
+#include <math.h>
+
+static void
+mapc2p_3x(double t, const double *xc, double* GKYL_RESTRICT xp, void *ctx)
+{
+  xp[0] = xc[0]; xp[1] = xc[1]; xp[2] = xc[2];
+}
+
+static void
+bfield_func_3x(double t, const double *xc, double* GKYL_RESTRICT fout, void *ctx)
+{
+  fout[0] = 0.0;
+  fout[1] = 0.0;
+  fout[2] = 1.0; // uniform Bz
+}
+
+// Shared fixture holding everything the moment constructor needs.
+struct fixture {
+  struct gkyl_basis basis, confBasis;
+  struct gkyl_range confLocal;
+  struct gk_geometry *gk_geom;
+  struct gkyl_velocity_map *gvm;
+  struct gkyl_position_map *pmap;
+};
+
+static void
+make_fixture(struct fixture *fx, int poly_order, int vdim)
+{
+  const int cdim = 1;
+  const int ndim = cdim + vdim;
+
+  double lower[GKYL_MAX_DIM], upper[GKYL_MAX_DIM];
+  int cells[GKYL_MAX_DIM];
+  lower[0] = -M_PI; upper[0] = M_PI; cells[0] = 4;
+  for (int d=0; d<vdim; d++) {
+    lower[cdim+d] = (d==0)? -2.0 : 0.0;
+    upper[cdim+d] = (d==0)?  2.0 : 2.0;
+    cells[cdim+d] = 2;
+  }
+
+  double confLower[1] = {lower[0]}, confUpper[1] = {upper[0]};
+  int confCells[1] = {cells[0]};
+  double velLower[3], velUpper[3];
+  int velCells[3];
+  for (int d=0; d<vdim; d++) {
+    velLower[d] = lower[cdim+d];
+    velUpper[d] = upper[cdim+d];
+    velCells[d] = cells[cdim+d];
+  }
+
+  struct gkyl_rect_grid grid, confGrid, velGrid;
+  gkyl_rect_grid_init(&grid, ndim, lower, upper, cells);
+  gkyl_rect_grid_init(&confGrid, cdim, confLower, confUpper, confCells);
+  gkyl_rect_grid_init(&velGrid, vdim, velLower, velUpper, velCells);
+
+  gkyl_cart_modal_serendip(&fx->confBasis, cdim, poly_order);
+  if (poly_order == 1)
+    gkyl_cart_modal_gkhybrid(&fx->basis, cdim, vdim);
+  else
+    gkyl_cart_modal_serendip(&fx->basis, ndim, poly_order);
+
+  int confGhost[] = {1, 1, 1};
+  struct gkyl_range confLocal_ext;
+  gkyl_create_grid_ranges(&confGrid, confGhost, &confLocal_ext, &fx->confLocal);
+
+  int velGhost[3] = {0};
+  struct gkyl_range velLocal, velLocal_ext;
+  gkyl_create_grid_ranges(&velGrid, velGhost, &velLocal_ext, &velLocal);
+  int ghost[GKYL_MAX_DIM] = {0};
+  ghost[0] = confGhost[0];
+  struct gkyl_range local, local_ext;
+  gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
+
+  fx->pmap = gkyl_position_map_null_new();
+
+  struct gkyl_gk_geometry_inp geometry_input = {
+    .geometry_id = GKYL_GEOMETRY_MAPC2P,
+    .world = {0.0, 0.0},
+    .mapc2p = mapc2p_3x,
+    .c2p_ctx = 0,
+    .bfield_func = bfield_func_3x,
+    .bfield_ctx = 0,
+    .position_map = fx->pmap,
+    .grid = confGrid,
+    .local = fx->confLocal,
+    .local_ext = confLocal_ext,
+    .global = fx->confLocal,
+    .global_ext = confLocal_ext,
+    .basis = fx->confBasis,
+  };
+  geometry_input.geo_grid = gkyl_gk_geometry_augment_grid(confGrid, geometry_input);
+  gkyl_create_grid_ranges(&geometry_input.geo_grid, confGhost,
+    &geometry_input.geo_local_ext, &geometry_input.geo_local);
+  gkyl_cart_modal_serendip(&geometry_input.geo_basis, 3, poly_order);
+  struct gk_geometry *gk_geom_3d = gkyl_gk_geometry_mapc2p_new(&geometry_input);
+  fx->gk_geom = gkyl_gk_geometry_deflate(gk_geom_3d, &geometry_input);
+  gkyl_gk_geometry_release(gk_geom_3d);
+
+  struct gkyl_mapc2p_inp c2p_in = { };
+  fx->gvm = gkyl_velocity_map_new(c2p_in, grid, velGrid, local, local_ext,
+    velLocal, velLocal_ext, false);
+}
+
+static void
+free_fixture(struct fixture *fx)
+{
+  gkyl_gk_geometry_release(fx->gk_geom);
+  gkyl_velocity_map_release(fx->gvm);
+  gkyl_position_map_release(fx->pmap);
+}
+
+static struct gkyl_mom_type*
+mk_mom(struct fixture *fx, enum gkyl_distribution_moments mt)
+{
+  return gkyl_mom_gyrokinetic_new(&fx->confBasis, &fx->basis, &fx->confLocal,
+    1.0, 1.0, fx->gvm, fx->gk_geom, NULL, mt, false);
+}
+
+void
+test_mom_types_1x2v()
+{
+  const int vdim = 2;
+  struct fixture fx;
+  make_fixture(&fx, 1, vdim);
+
+  struct { enum gkyl_distribution_moments mt; int nm; } cases[] = {
+    { GKYL_F_MOMENT_M0,             1 },
+    { GKYL_F_MOMENT_M1,             1 },
+    { GKYL_F_MOMENT_M2,             1 },
+    { GKYL_F_MOMENT_M2PAR,          1 },
+    { GKYL_F_MOMENT_M2PERP,         1 },
+    { GKYL_F_MOMENT_M3PAR,          1 },
+    { GKYL_F_MOMENT_M3PERP,         1 },
+    { GKYL_F_MOMENT_M0M1M2,         3 },
+    { GKYL_F_MOMENT_M0M1M2PARM2PERP, vdim+2 },
+    { GKYL_F_MOMENT_HAMILTONIAN,    3 },
+  };
+  int ncase = sizeof(cases)/sizeof(cases[0]);
+
+  for (int i=0; i<ncase; i++) {
+    struct gkyl_mom_type *m = mk_mom(&fx, cases[i].mt);
+    TEST_CHECK( m != NULL );
+    TEST_MSG("moment case index %d", i);
+    TEST_CHECK( m->cdim == 1 );
+    TEST_CHECK( m->pdim == 3 );
+    TEST_CHECK( m->poly_order == 1 );
+    TEST_CHECK( m->num_config == fx.confBasis.num_basis );
+    TEST_CHECK( m->num_phase == fx.basis.num_basis );
+    TEST_CHECK( m->num_mom == cases[i].nm );
+    TEST_CHECK( gkyl_mom_type_num_mom(m) == cases[i].nm );
+    gkyl_mom_type_release(m);
+  }
+
+  free_fixture(&fx);
+}
+
+void
+test_int_mom_types_1x2v()
+{
+  const int vdim = 2;
+  struct fixture fx;
+  make_fixture(&fx, 1, vdim);
+
+  // The integrated-moment constructor combines several moments into a single
+  // object; verify the multi-moment names produce the expected counts.
+  struct { enum gkyl_distribution_moments mt; int nm; } cases[] = {
+    { GKYL_F_MOMENT_M0M1M2,          3 },
+    { GKYL_F_MOMENT_M0M1M2PARM2PERP, vdim+2 },
+    { GKYL_F_MOMENT_HAMILTONIAN,     3 },
+  };
+  int ncase = sizeof(cases)/sizeof(cases[0]);
+
+  for (int i=0; i<ncase; i++) {
+    struct gkyl_mom_type *m = gkyl_int_mom_gyrokinetic_new(&fx.confBasis, &fx.basis,
+      &fx.confLocal, 1.0, 1.0, fx.gvm, fx.gk_geom, NULL, cases[i].mt, false);
+    TEST_CHECK( m != NULL );
+    TEST_MSG("int moment case index %d", i);
+    TEST_CHECK( m->cdim == 1 );
+    TEST_CHECK( m->pdim == 3 );
+    TEST_CHECK( m->num_mom == cases[i].nm );
+    gkyl_mom_type_release(m);
+  }
+
+  free_fixture(&fx);
+}
+
+void
+test_mom_types_1x1v()
+{
+  const int vdim = 1;
+  struct fixture fx;
+  make_fixture(&fx, 1, vdim);
+
+  // With a single velocity dimension, the par/perp split collapses; check a
+  // representative subset that is valid for 1v.
+  struct { enum gkyl_distribution_moments mt; int nm; } cases[] = {
+    { GKYL_F_MOMENT_M0,     1 },
+    { GKYL_F_MOMENT_M1,     1 },
+    { GKYL_F_MOMENT_M2,     1 },
+    { GKYL_F_MOMENT_M2PAR,  1 },
+    { GKYL_F_MOMENT_M3PAR,  1 },
+    { GKYL_F_MOMENT_M0M1M2, 3 },
+  };
+  int ncase = sizeof(cases)/sizeof(cases[0]);
+
+  for (int i=0; i<ncase; i++) {
+    struct gkyl_mom_type *m = mk_mom(&fx, cases[i].mt);
+    TEST_CHECK( m != NULL );
+    TEST_MSG("1v moment case index %d", i);
+    TEST_CHECK( m->cdim == 1 );
+    TEST_CHECK( m->pdim == 2 );
+    TEST_CHECK( m->num_mom == cases[i].nm );
+    gkyl_mom_type_release(m);
+  }
+
+  free_fixture(&fx);
+}
+
+TEST_LIST = {
+  { "mom_types_1x2v", test_mom_types_1x2v },
+  { "int_mom_types_1x2v", test_int_mom_types_1x2v },
+  { "mom_types_1x1v", test_mom_types_1x1v },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_nodal_ops_deflated.c b/gyrokinetic/unit/ctest_nodal_ops_deflated.c
new file mode 100644
index 0000000000..d96fd2a75c
--- /dev/null
+++ b/gyrokinetic/unit/ctest_nodal_ops_deflated.c
@@ -0,0 +1,306 @@
+#include <acutest.h>
+#include <math.h>
+#include <gkyl_array.h>
+#include <gkyl_array_rio.h>
+#include <gkyl_array_ops.h>
+#include <gkyl_eval_on_nodes.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_rect_decomp.h>
+#include <gkyl_util.h>
+#include <gkyl_basis.h>
+#include <gkyl_nodal_ops.h>
+#include <gkyl_deflate_zsurf.h>
+
+
+void check_same(struct gkyl_range range, struct gkyl_basis basis, struct gkyl_array *field1, struct gkyl_array* field2)
+{
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &range);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(&range, iter.idx);
+    const double *f1 = gkyl_array_cfetch(field1, lidx);
+    const double *f2 = gkyl_array_cfetch(field2, lidx);
+    for(int i = 0; i< basis.num_basis; i++)
+      TEST_CHECK( gkyl_compare(f1[i], f2[i], 1e-10) );
+  }
+}
+
+void
+proj_func(double t, const double *xn, double *fout, void *ctx)
+{
+  fout[0] = cos(xn[0])*sin(xn[1]);
+}
+
+void
+proj_func3d(double t, const double *xn, double *fout, void *ctx)
+{
+  fout[0] = cos(2*xn[0])*sin(xn[1])*xn[2]*xn[2]*xn[2];
+}
+
+void
+test_p1_deflated(){
+  // create  grid, ranges, basis
+  double lower[] = { 0.0, -1.5 }, upper[] = { 1.5, 1.5 };
+  int cells[] = { 8, 16 };
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 2, lower, upper, cells);
+  struct gkyl_range local, local_ext;
+  int nghost[GKYL_MAX_CDIM] = { 1, 1 };
+  gkyl_create_grid_ranges(&grid, nghost, &local_ext, &local);
+  int poly_order = 1;
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, 2, poly_order);
+  bool use_gpu = false;
+#ifdef GKYL_HAVE_CUDA
+  use_gpu = true;
+  struct gkyl_basis *basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
+  gkyl_cart_modal_serendip_cu_dev(basis_on_dev, 2, poly_order);
+#else
+  struct gkyl_basis *basis_on_dev = &basis;
+#endif
+
+  // Project Initial Function
+  struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+  gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func, 0);
+  gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+  gkyl_array_copy(funcdg_dev, funcdg);
+#else
+  struct gkyl_array *funcdg_dev = funcdg;
+#endif
+
+  // Construct Nodal range, field, and nodal ops object
+  int nodes[3] = { 1, 1, 1 };
+  for (int d=0; d<grid.ndim; ++d)
+    nodes[d] = grid.cells[d] + 1;
+  struct gkyl_range nrange;
+  gkyl_range_init_from_shape(&nrange, grid.ndim, nodes);
+  struct gkyl_array* nodal_fld = gkyl_array_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *nodal_fld_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
+#else
+  struct gkyl_array *nodal_fld_dev = nodal_fld;
+#endif
+  struct gkyl_nodal_ops *n2m = gkyl_nodal_ops_new(&basis, &grid, use_gpu);
+
+  // Create deflated 1d grid, ranges, basis, nodal range
+  struct gkyl_rect_grid deflated_grid;
+  struct gkyl_basis deflated_basis;
+  struct gkyl_range deflated_local, deflated_local_ext;
+  double deflated_lower[1] = { grid.lower[0]}, deflated_upper[1] = { grid.upper[0]};
+  int deflated_cells[1] = { grid.cells[0] };
+  gkyl_rect_grid_init(&deflated_grid, 1, deflated_lower, deflated_upper, deflated_cells);
+  int deflated_nghost[GKYL_MAX_CDIM] = { 1 };
+  gkyl_create_grid_ranges(&deflated_grid, deflated_nghost, &deflated_local_ext, &deflated_local);
+  gkyl_cart_modal_serendip(&deflated_basis, 1, poly_order);
+#ifdef GKYL_HAVE_CUDA
+  use_gpu = true;
+  struct gkyl_basis *deflated_basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
+  gkyl_cart_modal_serendip_cu_dev(deflated_basis_on_dev, deflated_grid.ndim, poly_order);
+#else
+  struct gkyl_basis *deflated_basis_on_dev = &deflated_basis;
+#endif
+  int deflated_nodes[1];
+  for (int d=0; d<deflated_grid.ndim; ++d)
+    deflated_nodes[d] = deflated_grid.cells[d] + 1;
+  struct gkyl_range deflated_nrange;
+  gkyl_range_init_from_shape(&deflated_nrange, deflated_grid.ndim, deflated_nodes);
+
+  // Create deflated modal field
+  struct gkyl_array *deflated_funcdg = gkyl_array_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *deflated_funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
+#else
+  struct gkyl_array *deflated_funcdg_dev = deflated_funcdg;
+#endif
+
+  // Created deflator objects and 1d nodal ops object
+  struct gkyl_deflate_zsurf *deflator_lo = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 0, use_gpu);
+  struct gkyl_deflate_zsurf *deflator_up = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 1, use_gpu);
+  struct gkyl_nodal_ops *n2m_1d = gkyl_nodal_ops_new(&deflated_basis, &deflated_grid, use_gpu);
+
+  // Loop along last dim (z) and call the 1d deflator to populate the 2d nodal field
+  for(int zidx = local.lower[1]; zidx <= local.upper[1]; zidx++){
+    gkyl_deflate_zsurf_advance(deflator_lo, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
+    gkyl_nodal_ops_m2n_deflated(n2m_1d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx-1);
+    if (zidx == local.upper[1]) {
+      gkyl_deflate_zsurf_advance(deflator_up, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
+      gkyl_nodal_ops_m2n_deflated(n2m_1d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx);
+    }
+  }
+
+  // Create output field
+  struct gkyl_array *funcdg2 = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *funcdg2_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+#else
+  struct gkyl_array *funcdg2_dev = funcdg2;
+#endif
+
+  // Transform back to modal
+  gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, false);
+#ifdef GKYL_HAVE_CUDA
+  gkyl_array_copy(funcdg2, funcdg2_dev);
+#endif
+  check_same(local, basis, funcdg, funcdg2);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
+
+#ifdef GKYL_HAVE_CUDA
+  gkyl_cu_free(basis_on_dev);
+  gkyl_array_release(funcdg_dev);
+  gkyl_array_release(deflated_funcdg_dev);
+  gkyl_array_release(funcdg2_dev);
+  gkyl_array_release(nodal_fld_dev);
+#endif
+  gkyl_array_release(funcdg);
+  gkyl_array_release(deflated_funcdg);
+  gkyl_array_release(funcdg2);
+  gkyl_array_release(nodal_fld);
+
+  gkyl_eval_on_nodes_release(eon);
+  gkyl_nodal_ops_release(n2m);
+  gkyl_nodal_ops_release(n2m_1d);
+  gkyl_deflate_zsurf_release(deflator_lo);
+  gkyl_deflate_zsurf_release(deflator_up);
+
+}
+
+void test_p1_deflated_3d(){
+  // create  grid, ranges, basis
+  double lower[] = { 0.0, -1.5, -2 }, upper[] = { 1.5, 1.5,2 };
+  int cells[] = { 8, 16 , 12};
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 3, lower, upper, cells);
+  struct gkyl_range local, local_ext;
+  int nghost[GKYL_MAX_CDIM] = { 1, 1 ,1};
+  gkyl_create_grid_ranges(&grid, nghost, &local_ext, &local);
+  int poly_order = 1;
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, 3, poly_order);
+  bool use_gpu = false;
+#ifdef GKYL_HAVE_CUDA
+  use_gpu = true;
+  struct gkyl_basis *basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
+  gkyl_cart_modal_serendip_cu_dev(basis_on_dev, 3, poly_order);
+#else
+  struct gkyl_basis *basis_on_dev = &basis;
+#endif
+
+  // Project Initial Function
+  struct gkyl_array *funcdg = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+  gkyl_eval_on_nodes *eon = gkyl_eval_on_nodes_new(&grid, &basis, 1, &proj_func3d, 0);
+  gkyl_eval_on_nodes_advance(eon, 0.0, &local, funcdg);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg, "proj_func.gkyl");
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+  gkyl_array_copy(funcdg_dev, funcdg);
+#else
+  struct gkyl_array *funcdg_dev = funcdg;
+#endif
+
+  // Construct Nodal range, field, and nodal ops object
+  int nodes[3] = { 1, 1, 1 };
+  for (int d=0; d<grid.ndim; ++d)
+    nodes[d] = grid.cells[d] + 1;
+  struct gkyl_range nrange;
+  gkyl_range_init_from_shape(&nrange, grid.ndim, nodes);
+  struct gkyl_array* nodal_fld = gkyl_array_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *nodal_fld_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, grid.ndim, nrange.volume);
+#else
+  struct gkyl_array *nodal_fld_dev = nodal_fld;
+#endif
+  struct gkyl_nodal_ops *n2m = gkyl_nodal_ops_new(&basis, &grid, use_gpu);
+
+  // Create deflated 2d grid, ranges, basis, nodal range
+  struct gkyl_rect_grid deflated_grid;
+  struct gkyl_basis deflated_basis;
+  struct gkyl_range deflated_local, deflated_local_ext;
+  double deflated_lower[2] = { grid.lower[0], grid.lower[1] };
+  double deflated_upper[2] = { grid.upper[0], grid.upper[1] };
+  int deflated_cells[2] = { grid.cells[0], grid.cells[1] };
+  gkyl_rect_grid_init(&deflated_grid, 2, deflated_lower, deflated_upper, deflated_cells);
+  int deflated_nghost[GKYL_MAX_CDIM] = { 1, 1 };
+  gkyl_create_grid_ranges(&deflated_grid, deflated_nghost, &deflated_local_ext, &deflated_local);
+  gkyl_cart_modal_serendip(&deflated_basis, 2, poly_order);
+#ifdef GKYL_HAVE_CUDA
+  use_gpu = true;
+  struct gkyl_basis *deflated_basis_on_dev = gkyl_cu_malloc(sizeof(struct gkyl_basis));
+  gkyl_cart_modal_serendip_cu_dev(deflated_basis_on_dev, deflated_grid.ndim, poly_order);
+#else
+  struct gkyl_basis *deflated_basis_on_dev = &deflated_basis;
+#endif
+  int deflated_nodes[2];
+  for (int d=0; d<deflated_grid.ndim; ++d)
+    deflated_nodes[d] = deflated_grid.cells[d] + 1;
+  struct gkyl_range deflated_nrange;
+  gkyl_range_init_from_shape(&deflated_nrange, deflated_grid.ndim, deflated_nodes);
+
+  // Create deflated modal field
+  struct gkyl_array *deflated_funcdg = gkyl_array_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *deflated_funcdg_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, deflated_basis.num_basis, local_ext.volume);
+#else
+  struct gkyl_array *deflated_funcdg_dev = deflated_funcdg;
+#endif
+
+  // Created deflator objects and 1d nodal ops object
+  struct gkyl_deflate_zsurf *deflator_lo = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 0, use_gpu);
+  struct gkyl_deflate_zsurf *deflator_up = gkyl_deflate_zsurf_new(&basis, &deflated_basis, 1, use_gpu);
+  struct gkyl_nodal_ops *n2m_2d = gkyl_nodal_ops_new(&deflated_basis, &deflated_grid, use_gpu);
+
+  // Loop along last dim (z) and call the 1d deflator to populate the 2d nodal field
+  for(int zidx = local.lower[2]; zidx <= local.upper[2]; zidx++){
+    gkyl_deflate_zsurf_advance(deflator_lo, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
+    gkyl_nodal_ops_m2n_deflated(n2m_2d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx-1);
+    if (zidx == local.upper[2]) {
+      gkyl_deflate_zsurf_advance(deflator_up, zidx, &local, &deflated_local, funcdg_dev, deflated_funcdg_dev, 1);
+      gkyl_nodal_ops_m2n_deflated(n2m_2d, deflated_basis_on_dev, &deflated_grid, &nrange, &deflated_nrange, &deflated_local, 1, nodal_fld_dev, deflated_funcdg_dev, zidx);
+    }
+  }
+
+  // Create output field
+  struct gkyl_array *funcdg2 = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+#ifdef GKYL_HAVE_CUDA
+  struct gkyl_array *funcdg2_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
+#else
+  struct gkyl_array *funcdg2_dev = funcdg2;
+#endif
+
+  // Transform back to modal
+  gkyl_nodal_ops_n2m(n2m, basis_on_dev, &grid, &nrange, &local, 1, nodal_fld_dev, funcdg2_dev, false);
+#ifdef GKYL_HAVE_CUDA
+  gkyl_array_copy(funcdg2, funcdg2_dev);
+#endif
+  check_same(local, basis, funcdg, funcdg2);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, funcdg2, "proj_func2.gkyl");
+
+#ifdef GKYL_HAVE_CUDA
+  gkyl_cu_free(basis_on_dev);
+  gkyl_array_release(funcdg_dev);
+  gkyl_array_release(deflated_funcdg_dev);
+  gkyl_array_release(funcdg2_dev);
+  gkyl_array_release(nodal_fld_dev);
+#endif
+  gkyl_array_release(funcdg);
+  gkyl_array_release(deflated_funcdg);
+  gkyl_array_release(funcdg2);
+  gkyl_array_release(nodal_fld);
+
+  gkyl_eval_on_nodes_release(eon);
+  gkyl_nodal_ops_release(n2m);
+  gkyl_nodal_ops_release(n2m_2d);
+  gkyl_deflate_zsurf_release(deflator_lo);
+  gkyl_deflate_zsurf_release(deflator_up);
+
+}
+
+
+TEST_LIST = {
+  { "test_p1_deflated", test_p1_deflated},
+  { "test_p1_deflated_3d", test_p1_deflated_3d},
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_position_map.c b/gyrokinetic/unit/ctest_position_map.c
index 7a12623471..d160d06008 100644
--- a/gyrokinetic/unit/ctest_position_map.c
+++ b/gyrokinetic/unit/ctest_position_map.c
@@ -4,6 +4,7 @@
 
 #include <gkyl_array.h>
 #include <gkyl_array_ops.h>
+#include <gkyl_basis.h>
 #include <gkyl_range.h>
 #include <gkyl_rect_decomp.h>
 #include <gkyl_rect_grid.h>
@@ -12,6 +13,7 @@
 #include <gkyl_position_map_priv.h>
 #include <gkyl_proj_on_basis.h>
 #include <gkyl_calc_bmag.h>
+#include <math.h>
 
 void
 test_nonuniform_position_map(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx)
@@ -572,16 +574,141 @@ test_position_map_numeric_calculate_1x()
   gkyl_array_release(bmag_global);
 }
 
+void
+test_null_map_is_identity()
+{
+  struct gkyl_position_map *pmap = gkyl_position_map_null_new();
+
+  TEST_CHECK( pmap != NULL );
+  TEST_CHECK( pmap->id == GKYL_PMAP_USER_INPUT );
+  TEST_CHECK( pmap->to_optimize == false );
+
+  // Each of the three maps must act as the identity, and each derivative as 1.
+  for (double z = -1.0; z <= 1.0; z += 0.25) {
+    for (int i=0; i<3; i++) {
+      double x[1] = {z}, y[1] = {0.0};
+      pmap->maps[i](0.0, x, y, pmap->ctxs[i]);
+      TEST_CHECK( gkyl_compare(y[0], z, 1e-15) );
+
+      double dy[1] = {0.0};
+      pmap->map_derivs[i](0.0, x, dy, pmap->ctxs[i]);
+      TEST_CHECK( gkyl_compare(dy[0], 1.0, 1e-15) );
+    }
+  }
+
+  gkyl_position_map_release(pmap);
+}
+
+static void
+nonuniform_map_1d(double t, const double *xn, double *fout, void *ctx)
+{
+  // A simple smooth monotone map z -> z + 0.1*sin(z) (identity-like near 0).
+  fout[0] = xn[0] + 0.1*sin(xn[0]);
+}
+
+void
+test_inew_constructor()
+{
+  int cells[] = {16};
+  int poly_order = 1;
+  double lower[] = {-1.0}, upper[] = {1.0};
+  int dim = 1;
+
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, dim, lower, upper, cells);
+  int ghost[] = {1};
+  struct gkyl_range local, local_ext;
+  gkyl_create_grid_ranges(&grid, ghost, &local_ext, &local);
+
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, dim, poly_order);
+
+  struct gkyl_position_map_inew_inp inp = {
+    .pmap_info = {
+      .maps = {nonuniform_map_1d, nonuniform_map_1d, nonuniform_map_1d},
+      .ctxs = {NULL, NULL, NULL},
+    },
+    .grid = grid,
+    .local = local, .local_ext = local_ext,
+    .global = local, .global_ext = local_ext,
+    .basis = basis,
+  };
+
+  struct gkyl_position_map *pmap = gkyl_position_map_inew(inp);
+
+  TEST_CHECK( pmap != NULL );
+  TEST_CHECK( pmap->grid.ndim == 1 );
+  TEST_CHECK( pmap->basis.poly_order == 1 );
+  TEST_CHECK( pmap->id == GKYL_PMAP_USER_INPUT );
+
+  // The user map must be wired in: check a couple of values.
+  double x[1] = {0.3}, y[1];
+  pmap->maps[0](0.0, x, y, pmap->ctxs[0]);
+  TEST_CHECK( gkyl_compare(y[0], 0.3 + 0.1*sin(0.3), 1e-14) );
+
+  gkyl_position_map_release(pmap);
+}
+
+void
+test_acquire_refcount()
+{
+  struct gkyl_position_map *pmap = gkyl_position_map_null_new();
+
+  struct gkyl_position_map *pmap2 = gkyl_position_map_acquire(pmap);
+  TEST_CHECK( pmap2 == pmap );
+
+  // Drop the first reference; the object must survive for use through pmap2.
+  gkyl_position_map_release(pmap);
+
+  double x[1] = {0.5}, y[1];
+  pmap2->maps[0](0.0, x, y, pmap2->ctxs[0]);
+  TEST_CHECK( gkyl_compare(y[0], 0.5, 1e-15) );
+
+  gkyl_position_map_release(pmap2);
+}
+
+void
+test_set_compression()
+{
+  struct gkyl_position_map *pmap = gkyl_position_map_null_new();
+
+  // With both compression factors zero the maps stay at the identity backups,
+  // but the geometric parameters must be recorded.
+  pmap->xpt_ctx->compression_factor = 0.0;
+  pmap->xpt_ctx->radial_compression_factor = 0.0;
+
+  double zcut = 1.25, zcenter = 0.1, w = 0.05, psisep = 0.7;
+  gkyl_position_map_set_compression(pmap, zcut, zcenter, w, psisep);
+
+  TEST_CHECK( gkyl_compare(pmap->xpt_ctx->zcut, zcut, 1e-15) );
+  TEST_CHECK( gkyl_compare(pmap->xpt_ctx->zcenter, zcenter, 1e-15) );
+  TEST_CHECK( gkyl_compare(pmap->xpt_ctx->w, w, 1e-15) );
+  TEST_CHECK( gkyl_compare(pmap->xpt_ctx->psisep, psisep, 1e-15) );
+
+  // Maps 0 and 1 fall back to identity backups when compression is disabled.
+  double x[1] = {0.4}, y[1];
+  pmap->maps[0](0.0, x, y, pmap->ctxs[0]);
+  TEST_CHECK( gkyl_compare(y[0], 0.4, 1e-15) );
+  pmap->maps[1](0.0, x, y, pmap->ctxs[1]);
+  TEST_CHECK( gkyl_compare(y[0], 0.4, 1e-15) );
+
+  gkyl_position_map_release(pmap);
+}
+
 TEST_LIST = {
   { "test_position_map_init_1x", test_position_map_init_1x },
   { "test_position_map_init_1x_null", test_position_map_init_1x_null },
   { "test_position_map_init_2x", test_position_map_init_2x },
   { "test_position_map_init_3x", test_position_map_init_3x },
   { "test_position_map_set", test_position_map_set },
-  { "test_gkyl_position_map_eval_mc2nu", test_gkyl_position_map_eval_mc2nu }, 
-  { "test_gkyl_position_map_slope", test_gkyl_position_map_slope },
-  { "test_position_polynomial_map_optimize_1x", test_position_polynomial_map_optimize_1x },
+  { "test_position_map_eval_mc2nu", test_gkyl_position_map_eval_mc2nu },
+  { "test_position_map_slope", test_gkyl_position_map_slope },
+  { "test_position_map_polynomial_optimize_1x", test_position_polynomial_map_optimize_1x },
   { "test_position_map_numeric_optimize_1x", test_position_map_numeric_optimize_1x },
   { "test_position_map_numeric_calculate_1x", test_position_map_numeric_calculate_1x },
+  { "test_position_map_null_map_is_identity", test_null_map_is_identity },
+  { "test_position_map_inew_constructor", test_inew_constructor },
+  { "test_position_map_acquire_refcount", test_acquire_refcount },
+  { "test_position_map_set_compression", test_set_compression },
   { NULL, NULL },
-};
\ No newline at end of file
+};
diff --git a/gyrokinetic/unit/ctest_prim_lbo_gyrokinetic.c b/gyrokinetic/unit/ctest_prim_lbo_gyrokinetic.c
new file mode 100644
index 0000000000..901b0ffc4d
--- /dev/null
+++ b/gyrokinetic/unit/ctest_prim_lbo_gyrokinetic.c
@@ -0,0 +1,80 @@
+// Test construction of the gyrokinetic LBO primitive-moment type object.
+// Verifies that the returned gkyl_prim_lbo_type has the expected dimension,
+// basis-count and udim fields for a variety of (cdim, vdim, poly_order)
+// configurations, and that acquire/release reference counting works.
+#include <acutest.h>
+
+#include <gkyl_basis.h>
+#include <gkyl_prim_lbo_type.h>
+#include <gkyl_prim_lbo_gyrokinetic.h>
+
+static void
+check_prim(int cdim, int vdim, int poly_order)
+{
+  int pdim = cdim + vdim;
+
+  struct gkyl_basis cbasis, pbasis;
+  gkyl_cart_modal_serendip(&cbasis, cdim, poly_order);
+  // Gyrokinetic phase space uses the gkhybrid basis at p=1, serendip otherwise.
+  if (poly_order == 1)
+    gkyl_cart_modal_gkhybrid(&pbasis, cdim, vdim);
+  else
+    gkyl_cart_modal_serendip(&pbasis, pdim, poly_order);
+
+  struct gkyl_prim_lbo_type *prim = gkyl_prim_lbo_gyrokinetic_new(&cbasis, &pbasis, false);
+
+  TEST_CHECK( prim != NULL );
+  TEST_CHECK( prim->cdim == cdim );
+  TEST_CHECK( prim->pdim == pdim );
+  TEST_CHECK( prim->poly_order == poly_order );
+  TEST_CHECK( prim->num_config == cbasis.num_basis );
+  TEST_CHECK( prim->num_phase == pbasis.num_basis );
+  // Gyrokinetic flow has a single (parallel) velocity component.
+  TEST_CHECK( prim->udim == 1 );
+  // Kernels must be wired up.
+  TEST_CHECK( prim->self_prim != NULL );
+  TEST_CHECK( prim->cross_prim != NULL );
+
+  gkyl_prim_lbo_type_release(prim);
+}
+
+void test_prim_ctor_1x1v_p1() { check_prim(1, 1, 1); }
+void test_prim_ctor_1x2v_p1() { check_prim(1, 2, 1); }
+void test_prim_ctor_2x2v_p1() { check_prim(2, 2, 1); }
+void test_prim_ctor_3x2v_p1() { check_prim(3, 2, 1); }
+void test_prim_ctor_1x1v_p2() { check_prim(1, 1, 2); }
+void test_prim_ctor_2x2v_p2() { check_prim(2, 2, 2); }
+
+void
+test_prim_ctor_acquire()
+{
+  struct gkyl_basis cbasis, pbasis;
+  gkyl_cart_modal_serendip(&cbasis, 1, 1);
+  gkyl_cart_modal_gkhybrid(&pbasis, 1, 2);
+
+  struct gkyl_prim_lbo_type *prim = gkyl_prim_lbo_gyrokinetic_new(&cbasis, &pbasis, false);
+
+  // Acquire a second reference; the underlying object must survive the first
+  // release, and report identical fields.
+  struct gkyl_prim_lbo_type *prim2 = gkyl_prim_lbo_type_acquire(prim);
+  TEST_CHECK( prim2 == prim );
+
+  gkyl_prim_lbo_type_release(prim);
+
+  TEST_CHECK( prim2->cdim == 1 );
+  TEST_CHECK( prim2->pdim == 3 );
+  TEST_CHECK( prim2->udim == 1 );
+
+  gkyl_prim_lbo_type_release(prim2);
+}
+
+TEST_LIST = {
+  { "prim_ctor_1x1v_p1", test_prim_ctor_1x1v_p1 },
+  { "prim_ctor_1x2v_p1", test_prim_ctor_1x2v_p1 },
+  { "prim_ctor_2x2v_p1", test_prim_ctor_2x2v_p1 },
+  { "prim_ctor_3x2v_p1", test_prim_ctor_3x2v_p1 },
+  { "prim_ctor_1x1v_p2", test_prim_ctor_1x1v_p2 },
+  { "prim_ctor_2x2v_p2", test_prim_ctor_2x2v_p2 },
+  { "prim_ctor_acquire", test_prim_ctor_acquire },
+  { NULL, NULL },
+};
diff --git a/gyrokinetic/unit/ctest_proj_gk_maxwellian_on_basis.c b/gyrokinetic/unit/ctest_proj_gk_maxwellian_on_basis.c
index 809b56b78a..15d94b6d18 100644
--- a/gyrokinetic/unit/ctest_proj_gk_maxwellian_on_basis.c
+++ b/gyrokinetic/unit/ctest_proj_gk_maxwellian_on_basis.c
@@ -263,7 +263,7 @@ test_1x2v_gk(int poly_order, bool use_gpu)
   else {
     sprintf(fname, "ctest_proj_gkmaxwellian_on_basis_prim_mom_1x2v_p%d_cpu.gkyl", poly_order);
   }   
-  gkyl_grid_sub_array_write(&grid, &local, 0, distf, fname);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, distf, fname);
 
   gkyl_array_release(den); 
   gkyl_array_release(udrift); 
@@ -551,7 +551,7 @@ test_3x2v_gk(int poly_order, bool use_gpu)
   else {
     sprintf(fname, "ctest_proj_gkmaxwellian_on_basis_prim_mom_3x2v_p%d_cpu.gkyl", poly_order);
   }   
-  gkyl_grid_sub_array_write(&grid, &local, 0, distf, fname);
+  // gkyl_grid_sub_array_write(&grid, &local, 0, distf, fname);
 
   // Calculate the moments and copy from device to host.
   struct gkyl_gk_maxwellian_moments_inp inp_calc = {
@@ -589,7 +589,7 @@ test_3x2v_gk(int poly_order, bool use_gpu)
   else {
     sprintf(fname_moms, "ctest_proj_gkmaxwellian_on_basis_prim_mom_3x2v_p%d_moms_cpu.gkyl", poly_order);
   }   
-  gkyl_grid_sub_array_write(&confGrid, &confLocal, 0, moms, fname_moms);
+  // gkyl_grid_sub_array_write(&confGrid, &confLocal, 0, moms, fname_moms);
 
   gkyl_array_release(den); 
   gkyl_array_release(udrift); 
diff --git a/gyrokinetic/unit/ctest_time_roots.c b/gyrokinetic/unit/ctest_time_roots.c
index 7c5a7d2809..a73863ee2e 100644
--- a/gyrokinetic/unit/ctest_time_roots.c
+++ b/gyrokinetic/unit/ctest_time_roots.c
@@ -216,8 +216,8 @@ compare_quad_and_cub(void)
   );
   struct gkyl_array *psi_cubic_DG = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, local_ext.volume);
   gkyl_proj_on_basis_advance(projCub, 0.0, &local, psi_cubic_DG);
-  gkyl_grid_sub_array_write(&grid, &local, 0, psi_cubic_DG, "psi_cubic.gkyl");
-  gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "psi_quad.gkyl");
+  // gkyl_grid_sub_array_write(&grid, &local, 0, psi_cubic_DG, "psi_cubic.gkyl");
+  // gkyl_grid_sub_array_write(&efit->rzgrid, &efit->rzlocal, 0, efit->psizr, "psi_quad.gkyl");
 
   // Now pick a value of Z. Let's choose Z = 0.0
   // We want to see how long each one takes to find the roots
diff --git a/gyrokinetic/unit/mctest_multib_allgather.c b/gyrokinetic/unit/mctest_multib_allgather.c
index 9771a1fd2e..7722de0e67 100644
--- a/gyrokinetic/unit/mctest_multib_allgather.c
+++ b/gyrokinetic/unit/mctest_multib_allgather.c
@@ -1037,7 +1037,7 @@ test_L_domain_allgather_dir0_cuts2_par()
     gkyl_rect_grid_init(&grid, 2, gridlo, gridup, cells);
     char str[50];
     sprintf(str, "lb%d_r%d.gkyl",bI, my_rank);
-    gkyl_grid_sub_array_write(&grid, global_ranges[bI], 0, array_global[bI], str);
+    // gkyl_grid_sub_array_write(&grid, global_ranges[bI], 0, array_global[bI], str);
   }
 
   for (int bI = 0; bI<num_local_blocks; ++bI) {
@@ -1277,7 +1277,7 @@ test_SOL_domain_allgather_dir1_cuts2_par(bool use_gpu)
     gkyl_rect_grid_init(&grid, 2, gridlo, gridup, cells);
     char str[50];
     sprintf(str, "lb%d_r%d.gkyl",bI, my_rank);
-    gkyl_grid_sub_array_write(&grid, global_ranges[bI], 0, array_global_ho[bI], str);
+    // gkyl_grid_sub_array_write(&grid, global_ranges[bI], 0, array_global_ho[bI], str);
   }
 
   printf("checking\n");
diff --git a/gyrokinetic/zero/gk_geometry.c b/gyrokinetic/zero/gk_geometry.c
index 1b944537e3..ae3ed65a20 100644
--- a/gyrokinetic/zero/gk_geometry.c
+++ b/gyrokinetic/zero/gk_geometry.c
@@ -70,7 +70,6 @@ gkyl_gk_geometry_new(struct gk_geometry* geo_host, struct gkyl_gk_geometry_inp *
     gkyl_cart_modal_serendip(&up->surf_basis, 1, 0);
   }
   up->num_surf_basis = up->surf_basis.num_basis;
-  up->num_surf_basis = up->surf_basis.num_basis;
 
   gk_geometry_corn_alloc_expansions(up);
   gk_geometry_corn_alloc_nodal(up);