diff --git a/core/unit/ctest_proj_on_basis.c b/core/unit/ctest_proj_on_basis.c
index ae780a95a..4ba855058 100644
--- a/core/unit/ctest_proj_on_basis.c
+++ b/core/unit/ctest_proj_on_basis.c
@@ -1,9 +1,11 @@
 #include <acutest.h>
 
+#include <gkyl_alloc.h>
 #include <gkyl_proj_on_basis.h>
 #include <gkyl_range.h>
 #include <gkyl_rect_decomp.h>
 #include <gkyl_rect_grid.h>
+#include <gkyl_util.h>
 #include <math.h>
 
 void evalFunc(double t, const double *xn, double* restrict fout, void *ctx)
@@ -380,11 +382,221 @@ test_3_3d()
   gkyl_array_release(distf);
 }
 
+// Cuda specific tests
+#ifdef GKYL_HAVE_CUDA
+
+// Getters for device function addresses, defined in ctest_proj_on_basis_cu.cu.
+// The device functions must match their host counterparts in this file.
+evalf_t ctest_proj_on_basis_f_1d_cu_dev_ptr(void);
+evalf_t ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void);
+proj_on_basis_c2p_t ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void);
+
+// Context for the 2d two-component function; must match the definition in
+// ctest_proj_on_basis_cu.cu.
+struct ctest_proj_on_basis_2d_ctx {
+  double c0, c1;
+};
+
+// Host counterpart of the device function ctest_pob_f_2d_2c.
+void evalFunc_2d_2c(double t, const double *xn, double* restrict fout, void *ctx)
+{
+  struct ctest_proj_on_basis_2d_ctx *tctx = ctx;
+  double x = xn[0], y = xn[1];
+  fout[0] = tctx->c0 + x*y;
+  fout[1] = tctx->c1*x*x + y;
+}
+
+// Host counterpart of the device function ctest_pob_c2p_1d.
+void c2pFunc_1d(const double *xcomp, double *xphys, void *ctx)
+{
+  xphys[0] = 0.5*xcomp[0] + 0.1;
+}
+
+// Compare a device array against a host reference, coefficient by coefficient.
+static void
+check_same_dev_ho(const struct gkyl_array *ref_ho, struct gkyl_array *arr_cu)
+{
+  struct gkyl_array *arr_ho = gkyl_array_new(GKYL_DOUBLE, arr_cu->ncomp, arr_cu->size);
+  gkyl_array_copy(arr_ho, arr_cu);
+
+  for (size_t i=0; i<ref_ho->size; ++i) {
+    const double *ref_c = gkyl_array_cfetch(ref_ho, i);
+    const double *arr_c = gkyl_array_cfetch(arr_ho, i);
+    for (size_t k=0; k<ref_ho->ncomp; ++k) {
+      TEST_CHECK( gkyl_compare(ref_c[k], arr_c[k], 1e-14) );
+      TEST_MSG("cell %zu coeff %zu | Expected: %.13e | Produced: %.13e", i, k, ref_c[k], arr_c[k]);
+    }
+  }
+
+  gkyl_array_release(arr_ho);
+}
+
+void
+test_1_cu()
+{
+  int poly_order = 1;
+  double lower[] = {-2.0}, upper[] = {2.0};
+  int cells[] = {2};
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, 1, poly_order);
+
+  int nghost[GKYL_MAX_DIM] = { 0 };
+  struct gkyl_range arr_range, arr_ext_range;
+  gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range);
+
+  // Reference: project on the host.
+  gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis,
+    poly_order+1, 1, evalFunc, NULL);
+  struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf);
+
+  // Project the same function on the device.
+  gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+      .grid = &grid,
+      .basis = &basis,
+      .qtype = GKYL_GAUSS_QUAD,
+      .num_quad = poly_order+1,
+      .num_ret_vals = 1,
+      .eval = ctest_proj_on_basis_f_1d_cu_dev_ptr(),
+      .ctx = NULL,
+      .use_gpu = true,
+    }
+  );
+  struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu);
+
+  check_same_dev_ho(distf, distf_cu);
+
+  gkyl_proj_on_basis_release(projDistf);
+  gkyl_proj_on_basis_release(projDistf_cu);
+  gkyl_array_release(distf);
+  gkyl_array_release(distf_cu);
+}
+
+void
+test_2d_2c_cu()
+{
+  int poly_order = 2;
+  double lower[] = {-2.0,-1.0}, upper[] = {2.0,3.0};
+  int cells[] = {4, 3};
+  int ndim = sizeof(cells)/sizeof(cells[0]);
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, ndim, lower, upper, cells);
+
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, ndim, poly_order);
+
+  int nghost[GKYL_MAX_DIM] = { 0 };
+  struct gkyl_range arr_range, arr_ext_range;
+  gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range);
+
+  int num_ret_vals = 2;
+  struct ctest_proj_on_basis_2d_ctx tctx = { .c0 = 0.5, .c1 = 2.0 };
+
+  // Reference: project on the host.
+  gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis,
+    poly_order+1, num_ret_vals, evalFunc_2d_2c, &tctx);
+  struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE,
+    num_ret_vals*basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf);
+
+  // Project the same function on the device; the context must be placed
+  // in device memory by the user.
+  struct ctest_proj_on_basis_2d_ctx *tctx_cu = gkyl_cu_malloc(sizeof(struct ctest_proj_on_basis_2d_ctx));
+  gkyl_cu_memcpy(tctx_cu, &tctx, sizeof(struct ctest_proj_on_basis_2d_ctx), GKYL_CU_MEMCPY_H2D);
+
+  gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+      .grid = &grid,
+      .basis = &basis,
+      .qtype = GKYL_GAUSS_QUAD,
+      .num_quad = poly_order+1,
+      .num_ret_vals = num_ret_vals,
+      .eval = ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(),
+      .ctx = tctx_cu,
+      .use_gpu = true,
+    }
+  );
+  struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE,
+    num_ret_vals*basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu);
+
+  check_same_dev_ho(distf, distf_cu);
+
+  gkyl_proj_on_basis_release(projDistf);
+  gkyl_proj_on_basis_release(projDistf_cu);
+  gkyl_array_release(distf);
+  gkyl_array_release(distf_cu);
+  gkyl_cu_free(tctx_cu);
+}
+
+void
+test_c2p_1d_cu()
+{
+  int poly_order = 1;
+  double lower[] = {-2.0}, upper[] = {2.0};
+  int cells[] = {8};
+  struct gkyl_rect_grid grid;
+  gkyl_rect_grid_init(&grid, 1, lower, upper, cells);
+
+  struct gkyl_basis basis;
+  gkyl_cart_modal_serendip(&basis, 1, poly_order);
+
+  int nghost[GKYL_MAX_DIM] = { 0 };
+  struct gkyl_range arr_range, arr_ext_range;
+  gkyl_create_grid_ranges(&grid, nghost, &arr_ext_range, &arr_range);
+
+  // Reference: project on the host with the host c2p mapping.
+  gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+      .grid = &grid,
+      .basis = &basis,
+      .qtype = GKYL_GAUSS_QUAD,
+      .num_quad = poly_order+1,
+      .num_ret_vals = 1,
+      .eval = evalFunc,
+      .c2p_func = c2pFunc_1d,
+    }
+  );
+  struct gkyl_array *distf = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf, 0.0, &arr_range, distf);
+
+  // Project on the device with the equivalent device c2p mapping.
+  gkyl_proj_on_basis *projDistf_cu = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+      .grid = &grid,
+      .basis = &basis,
+      .qtype = GKYL_GAUSS_QUAD,
+      .num_quad = poly_order+1,
+      .num_ret_vals = 1,
+      .eval = ctest_proj_on_basis_f_1d_cu_dev_ptr(),
+      .c2p_func = ctest_proj_on_basis_c2p_1d_cu_dev_ptr(),
+      .use_gpu = true,
+    }
+  );
+  struct gkyl_array *distf_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, basis.num_basis, arr_range.volume);
+  gkyl_proj_on_basis_advance(projDistf_cu, 0.0, &arr_range, distf_cu);
+
+  check_same_dev_ho(distf, distf_cu);
+
+  gkyl_proj_on_basis_release(projDistf);
+  gkyl_proj_on_basis_release(projDistf_cu);
+  gkyl_array_release(distf);
+  gkyl_array_release(distf_cu);
+}
+
+#endif
+
 TEST_LIST = {
   { "test_1", test_1 },
-  { "test_2", test_2 },  
-  { "test_2_2d", test_2_2d },  
-  { "test_2_3d", test_2_3d },  
-  { "test_3_3d", test_3_3d },  
+  { "test_2", test_2 },
+  { "test_2_2d", test_2_2d },
+  { "test_2_3d", test_2_3d },
+  { "test_3_3d", test_3_3d },
+#ifdef GKYL_HAVE_CUDA
+  { "test_1_cu", test_1_cu },
+  { "test_2d_2c_cu", test_2d_2c_cu },
+  { "test_c2p_1d_cu", test_c2p_1d_cu },
+#endif
   { NULL, NULL },
 };
diff --git a/core/unit/ctest_proj_on_basis_cu.cu b/core/unit/ctest_proj_on_basis_cu.cu
new file mode 100644
index 000000000..476bad959
--- /dev/null
+++ b/core/unit/ctest_proj_on_basis_cu.cu
@@ -0,0 +1,64 @@
+/* -*- c++ -*- */
+
+// Device-side functions (and getters for their device addresses) used by the
+// GPU tests in ctest_proj_on_basis.c. Each function must match its host
+// counterpart in that file exactly.
+
+extern "C" {
+#include <gkyl_proj_on_basis.h>
+#include <gkyl_util.h>
+
+evalf_t ctest_proj_on_basis_f_1d_cu_dev_ptr(void);
+evalf_t ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void);
+proj_on_basis_c2p_t ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void);
+}
+
+// Context for the 2d two-component function; must match the definition in
+// ctest_proj_on_basis.c.
+struct ctest_proj_on_basis_2d_ctx {
+  double c0, c1;
+};
+
+GKYL_CU_DH static void
+ctest_pob_f_1d(double t, const double *xn, double *fout, void *ctx)
+{
+  double x = xn[0];
+  fout[0] = x*x;
+}
+
+GKYL_CU_DH static void
+ctest_pob_f_2d_2c(double t, const double *xn, double *fout, void *ctx)
+{
+  struct ctest_proj_on_basis_2d_ctx *tctx = (struct ctest_proj_on_basis_2d_ctx *) ctx;
+  double x = xn[0], y = xn[1];
+  fout[0] = tctx->c0 + x*y;
+  fout[1] = tctx->c1*x*x + y;
+}
+
+GKYL_CU_DH static void
+ctest_pob_c2p_1d(const double *xcomp, double *xphys, void *ctx)
+{
+  xphys[0] = 0.5*xcomp[0] + 0.1;
+}
+
+GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_f_1d, evalf_t, ctest_pob_f_1d_getter);
+GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_f_2d_2c, evalf_t, ctest_pob_f_2d_2c_getter);
+GKYL_DEFINE_CU_DEV_FUNC_GETTER(ctest_pob_c2p_1d, proj_on_basis_c2p_t, ctest_pob_c2p_1d_getter);
+
+extern "C" evalf_t
+ctest_proj_on_basis_f_1d_cu_dev_ptr(void)
+{
+  return ctest_pob_f_1d_getter();
+}
+
+extern "C" evalf_t
+ctest_proj_on_basis_f_2d_2c_cu_dev_ptr(void)
+{
+  return ctest_pob_f_2d_2c_getter();
+}
+
+extern "C" proj_on_basis_c2p_t
+ctest_proj_on_basis_c2p_1d_cu_dev_ptr(void)
+{
+  return ctest_pob_c2p_1d_getter();
+}
diff --git a/core/zero/gkyl_proj_on_basis.h b/core/zero/gkyl_proj_on_basis.h
index 878bc6efe..993995ed4 100644
--- a/core/zero/gkyl_proj_on_basis.h
+++ b/core/zero/gkyl_proj_on_basis.h
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <stdbool.h>
+
 #include <gkyl_array.h>
 #include <gkyl_basis.h>
 #include <gkyl_eqn_type.h>
@@ -28,6 +30,13 @@ struct gkyl_proj_on_basis_inp {
   proj_on_basis_c2p_t c2p_func; // Function that transforms a set of ndim
                                 // computational coordinates to physical ones.
   void *c2p_func_ctx; // Context for c2p_func.
+
+  bool use_gpu; // Whether to run the projection on the GPU. If true, 'eval'
+                // (and 'c2p_func', if provided) must be device function
+                // pointers (see GKYL_DEFINE_CU_DEV_FUNC_GETTER in gkyl_util.h)
+                // and 'ctx'/'c2p_func_ctx' must point to GPU-resident memory;
+                // the user must place their context on the GPU themselves as
+                // this updater has no way to perform that copy.
 };
 
 /**
@@ -63,7 +72,7 @@ gkyl_proj_on_basis *gkyl_proj_on_basis_new(const struct gkyl_rect_grid *grid,
  * @param pob Project on basis updater to run
  * @param tm Time at which projection must be computed
  * @param update_rng Range on which to run projection.
- * @param out Output array
+ * @param out Output array (a device array if the updater was created with use_gpu=true).
  */
 void gkyl_proj_on_basis_advance(const gkyl_proj_on_basis *pob,
   double tm, const struct gkyl_range *update_rng, struct gkyl_array *out);
diff --git a/core/zero/gkyl_proj_on_basis_priv.h b/core/zero/gkyl_proj_on_basis_priv.h
new file mode 100644
index 000000000..5bedc7a37
--- /dev/null
+++ b/core/zero/gkyl_proj_on_basis_priv.h
@@ -0,0 +1,68 @@
+// Private header for the proj_on_basis updater. Not for direct use in user code.
+#pragma once
+
+#include <gkyl_array.h>
+#include <gkyl_proj_on_basis.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_util.h>
+
+struct gkyl_proj_on_basis {
+  struct gkyl_rect_grid grid;
+  int num_quad; // number of quadrature points to use in each direction
+  int num_ret_vals; // number of values returned by eval function
+  evalf_t eval; // function to project
+  void *ctx; // evaluation context
+
+  int num_basis; // number of basis functions
+  int tot_quad; // total number of quadrature points
+  struct gkyl_array *ordinates; // ordinates for quadrature
+  struct gkyl_array *weights; // weights for quadrature
+  struct gkyl_array *basis_at_ords; // basis functions at ordinates
+
+  proj_on_basis_c2p_t c2p; // Function transformin comp to phys coords.
+  void *c2p_ctx; // Context for the c2p mapping.
+
+  bool use_gpu; // Whether the projection runs on the GPU.
+  // Device copies of the quadrature data (only allocated when use_gpu=true).
+  struct gkyl_array *ordinates_cu;
+  struct gkyl_array *weights_cu;
+  struct gkyl_array *basis_at_ords_cu;
+  struct gkyl_proj_on_basis *on_dev; // Device clone of this updater (points to itself on CPU).
+};
+
+GKYL_CU_DH
+static inline void
+proj_on_basis_log_to_comp(int ndim, const double *eta,
+  const double * GKYL_RESTRICT dx, const double * GKYL_RESTRICT xc,
+  double* GKYL_RESTRICT xout)
+{
+  // Convert logical to computational coordinates.
+  for (int d=0; d<ndim; ++d) xout[d] = 0.5*dx[d]*eta[d]+xc[d];
+}
+
+#ifdef GKYL_HAVE_CUDA
+
+/**
+ * Create the device clone of the updater. The host-side updater must already
+ * have its quadrature device mirrors (ordinates_cu, weights_cu,
+ * basis_at_ords_cu) allocated and populated, and its eval/c2p members must be
+ * device function pointers (c2p may be NULL, selecting an identity mapping
+ * assigned on the device).
+ *
+ * @param up Host-side updater.
+ * @return Pointer to the device clone.
+ */
+struct gkyl_proj_on_basis* gkyl_proj_on_basis_cu_dev_new(struct gkyl_proj_on_basis *up);
+
+/**
+ * Run the projection on the GPU.
+ *
+ * @param up Project on basis updater (created with use_gpu=true).
+ * @param tm Time at which projection must be computed.
+ * @param update_range Range on which to run projection.
+ * @param arr Output array (must be a device array).
+ */
+void gkyl_proj_on_basis_advance_cu(const struct gkyl_proj_on_basis *up,
+  double tm, const struct gkyl_range *update_range, struct gkyl_array *arr);
+
+#endif
diff --git a/core/zero/gkyl_util.h b/core/zero/gkyl_util.h
index 8c434c122..83a437484 100644
--- a/core/zero/gkyl_util.h
+++ b/core/zero/gkyl_util.h
@@ -127,6 +127,37 @@ inline cudaError_t __checkCudaErrors__(cudaError_t code, const char *func, const
   return code;
 }
 
+#ifdef __CUDACC__
+// Define a host-callable function 'getter' that returns the device-side
+// address of the __device__ (or __device__ __host__) function 'func', so it
+// can be handed to updaters that call it inside kernels (e.g.
+// gkyl_proj_on_basis with use_gpu=true). The pointer is assigned inside a
+// kernel, where 'func' resolves to its device address (this avoids
+// troublesome cudaMemcpyFromSymbol). Only available in CUDA translation
+// units (.cu files); expose the getter to C code with an extern "C" wrapper.
+// Example:
+//
+//   GKYL_CU_DH static void ef(double t, const double *xn, double *fout, void *ctx) { ... }
+//   GKYL_DEFINE_CU_DEV_FUNC_GETTER(ef, evalf_t, ef_dev_ptr);
+//   ...
+//   evalf_t eval_dev = ef_dev_ptr();
+#define GKYL_DEFINE_CU_DEV_FUNC_GETTER(func, func_type, getter)         \
+  static __global__ void getter##_set_cu_ker(func_type *fptr_d)         \
+  {                                                                     \
+    *fptr_d = func;                                                     \
+  }                                                                     \
+  static func_type getter(void)                                         \
+  {                                                                     \
+    func_type *fptr_d, fptr_ho;                                         \
+    checkCuda(cudaMalloc(&fptr_d, sizeof(func_type)));                  \
+    getter##_set_cu_ker<<<1,1>>>(fptr_d);                               \
+    checkCuda(cudaMemcpy(&fptr_ho, fptr_d, sizeof(func_type),           \
+      cudaMemcpyDeviceToHost));                                         \
+    checkCuda(cudaFree(fptr_d));                                        \
+    return fptr_ho;                                                     \
+  }
+#endif // __CUDACC__
+
 #else
 
 #undef GKYL_HAVE_CUDA
diff --git a/core/zero/proj_on_basis.c b/core/zero/proj_on_basis.c
index 2d2303817..2f9bcb393 100644
--- a/core/zero/proj_on_basis.c
+++ b/core/zero/proj_on_basis.c
@@ -5,25 +5,9 @@
 #include <gkyl_array.h>
 #include <gkyl_gauss_quad_data.h>
 #include <gkyl_proj_on_basis.h>
+#include <gkyl_proj_on_basis_priv.h>
 #include <gkyl_range.h>
 
-struct gkyl_proj_on_basis {
-  struct gkyl_rect_grid grid;
-  int num_quad; // number of quadrature points to use in each direction
-  int num_ret_vals; // number of values returned by eval function
-  evalf_t eval; // function to project
-  void *ctx; // evaluation context
-
-  int num_basis; // number of basis functions
-  int tot_quad; // total number of quadrature points
-  struct gkyl_array *ordinates; // ordinates for quadrature
-  struct gkyl_array *weights; // weights for quadrature
-  struct gkyl_array *basis_at_ords; // basis functions at ordinates
-
-  proj_on_basis_c2p_t c2p; // Function transformin comp to phys coords.
-  void *c2p_ctx; // Context for the c2p mapping.
-};
-
 // Identity comp to phys coord mapping, for when user doesn't provide a map.
 static inline void
 c2p_identity(const double *xcomp, double *xphys, void *ctx)
@@ -47,6 +31,7 @@ gkyl_proj_on_basis_new(const struct gkyl_rect_grid *grid, const struct gkyl_basi
       .ctx = ctx,
       .c2p_func = 0,
       .c2p_func_ctx = NULL,
+      .use_gpu = false,
     }
   );
 }
@@ -62,8 +47,15 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp)
   up->eval = inp->eval;
   up->ctx = inp->ctx;
   up->num_basis = inp->basis->num_basis;
+  up->use_gpu = inp->use_gpu;
 
-  if (inp->c2p_func == 0) {
+  if (up->use_gpu) {
+    // On the GPU eval and c2p_func must be device function pointers; a NULL
+    // c2p_func selects an identity mapping assigned on the device.
+    up->c2p = inp->c2p_func;
+    up->c2p_ctx = inp->c2p_func_ctx;
+  }
+  else if (inp->c2p_func == 0) {
     up->c2p = c2p_identity;
     up->c2p_ctx = &up->grid; // Use grid as the context since all we need is ndim.
   }
@@ -87,7 +79,7 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp)
   }
   else if (inp->qtype == GKYL_GAUSS_LOBATTO_QUAD) {
     assert( (num_quad > 1) && (num_quad <= gkyl_gauss_max) );
-    
+
     // Gauss-Lobatto quadrature
     memcpy(ordinates1, gkyl_gauss_lobatto_ordinates[num_quad], sizeof(double[num_quad]));
     memcpy(weights1, gkyl_gauss_lobatto_weights[num_quad], sizeof(double[num_quad]));
@@ -105,7 +97,7 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp)
 
   int tot_quad = up->tot_quad = qrange.volume;
 
-  // create ordinates and weights for multi-D quadrature  
+  // create ordinates and weights for multi-D quadrature
   up->ordinates = gkyl_array_new(GKYL_DOUBLE, inp->grid->ndim, tot_quad);
   up->weights = gkyl_array_new(GKYL_DOUBLE, 1, tot_quad);
 
@@ -114,12 +106,12 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp)
 
   while (gkyl_range_iter_next(&iter)) {
     long node = gkyl_range_idx(&qrange, iter.idx);
-    
+
     // set ordinates
     double *ord = gkyl_array_fetch(up->ordinates, node);
     for (int i=0; i<inp->grid->ndim; ++i)
       ord[i] = ordinates1[iter.idx[i]-qrange.lower[i]];
-    
+
     // set weights
     double *wgt = gkyl_array_fetch(up->weights, node);
     wgt[0] = 1.0;
@@ -133,6 +125,25 @@ gkyl_proj_on_basis_inew(const struct gkyl_proj_on_basis_inp *inp)
     inp->basis->eval(gkyl_array_fetch(up->ordinates, n),
       gkyl_array_fetch(up->basis_at_ords, n));
 
+  up->ordinates_cu = 0;
+  up->weights_cu = 0;
+  up->basis_at_ords_cu = 0;
+  up->on_dev = up; // On the CPU the updater points to itself.
+#ifdef GKYL_HAVE_CUDA
+  if (up->use_gpu) {
+    // Mirror the quadrature data on the device and create the device clone.
+    up->ordinates_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, inp->grid->ndim, tot_quad);
+    up->weights_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, 1, tot_quad);
+    up->basis_at_ords_cu = gkyl_array_cu_dev_new(GKYL_DOUBLE, inp->basis->num_basis, tot_quad);
+    gkyl_array_copy(up->ordinates_cu, up->ordinates);
+    gkyl_array_copy(up->weights_cu, up->weights);
+    gkyl_array_copy(up->basis_at_ords_cu, up->basis_at_ords);
+    up->on_dev = gkyl_proj_on_basis_cu_dev_new(up);
+  }
+#else
+  assert(!up->use_gpu);
+#endif
+
   return up;
 }
 
@@ -146,15 +157,6 @@ double* gkyl_proj_on_basis_fetch_ordinate(const struct gkyl_proj_on_basis *up, l
   return gkyl_array_fetch(up->ordinates, node);
 }
 
-static inline void
-log_to_comp(int ndim, const double *eta,
-  const double * GKYL_RESTRICT dx, const double * GKYL_RESTRICT xc,
-  double* GKYL_RESTRICT xout)
-{
-  // Convert logical to computational coordinates.
-  for (int d=0; d<ndim; ++d) xout[d] = 0.5*dx[d]*eta[d]+xc[d];
-}
-
 void
 gkyl_proj_on_basis_quad(const struct gkyl_proj_on_basis *up, const struct gkyl_array *fun_at_ords, double* f)
 {
@@ -186,20 +188,27 @@ void
 gkyl_proj_on_basis_advance(const struct gkyl_proj_on_basis *up,
   double tm, const struct gkyl_range *update_range, struct gkyl_array *arr)
 {
+#ifdef GKYL_HAVE_CUDA
+  if (up->use_gpu) {
+    gkyl_proj_on_basis_advance_cu(up, tm, update_range, arr);
+    return;
+  }
+#endif
+
   double xc[GKYL_MAX_DIM], xmu[GKYL_MAX_DIM];
 
   int num_ret_vals = up->num_ret_vals;
   int tot_quad = up->tot_quad;
   struct gkyl_array *fun_at_ords = gkyl_array_new(GKYL_DOUBLE, num_ret_vals, tot_quad);
-  
+
   struct gkyl_range_iter iter;
   gkyl_range_iter_init(&iter, update_range);
-  
+
   while (gkyl_range_iter_next(&iter)) {
     gkyl_rect_grid_cell_center(&up->grid, iter.idx, xc);
 
     for (int i=0; i<tot_quad; ++i) {
-      log_to_comp(up->grid.ndim, gkyl_array_cfetch(up->ordinates, i),
+      proj_on_basis_log_to_comp(up->grid.ndim, gkyl_array_cfetch(up->ordinates, i),
         up->grid.dx, xc, xmu);
       up->c2p(xmu, xmu, up->c2p_ctx);
       up->eval(tm, xmu, gkyl_array_fetch(fun_at_ords, i), up->ctx);
@@ -218,5 +227,13 @@ gkyl_proj_on_basis_release(struct gkyl_proj_on_basis* up)
   gkyl_array_release(up->ordinates);
   gkyl_array_release(up->weights);
   gkyl_array_release(up->basis_at_ords);
+#ifdef GKYL_HAVE_CUDA
+  if (up->use_gpu) {
+    gkyl_array_release(up->ordinates_cu);
+    gkyl_array_release(up->weights_cu);
+    gkyl_array_release(up->basis_at_ords_cu);
+    gkyl_cu_free(up->on_dev);
+  }
+#endif
   gkyl_free(up);
 }
diff --git a/core/zero/proj_on_basis_cu.cu b/core/zero/proj_on_basis_cu.cu
new file mode 100644
index 000000000..31b0fcc6b
--- /dev/null
+++ b/core/zero/proj_on_basis_cu.cu
@@ -0,0 +1,123 @@
+/* -*- c++ -*- */
+
+extern "C" {
+#include <assert.h>
+
+#include <gkyl_alloc.h>
+#include <gkyl_array.h>
+#include <gkyl_proj_on_basis.h>
+#include <gkyl_proj_on_basis_priv.h>
+#include <gkyl_range.h>
+#include <gkyl_rect_grid.h>
+#include <gkyl_util.h>
+}
+
+// Identity comp to phys coord mapping, for when user doesn't provide a map.
+GKYL_CU_D static void
+proj_on_basis_c2p_identity_cu(const double *xcomp, double *xphys, void *ctx)
+{
+  struct gkyl_rect_grid *grid = (struct gkyl_rect_grid *) ctx;
+  int ndim = grid->ndim;
+  for (int d=0; d<ndim; d++) xphys[d] = xcomp[d];
+}
+
+__global__ static void
+proj_on_basis_set_c2p_identity_cu_ker(struct gkyl_proj_on_basis *up)
+{
+  // Assigned in device code so the function pointer is a device address.
+  up->c2p = proj_on_basis_c2p_identity_cu;
+  up->c2p_ctx = &up->grid; // Use grid as the context since all we need is ndim.
+}
+
+struct gkyl_proj_on_basis*
+gkyl_proj_on_basis_cu_dev_new(struct gkyl_proj_on_basis *up)
+{
+  struct gkyl_proj_on_basis *up_cu = (struct gkyl_proj_on_basis*)
+    gkyl_cu_malloc(sizeof(struct gkyl_proj_on_basis));
+
+  // Clone with device pointers to the quadrature data. The eval/c2p members
+  // and their contexts are device pointers provided by the user.
+  struct gkyl_proj_on_basis up_ho = *up;
+  up_ho.ordinates = up->ordinates_cu->on_dev;
+  up_ho.weights = up->weights_cu->on_dev;
+  up_ho.basis_at_ords = up->basis_at_ords_cu->on_dev;
+  up_ho.on_dev = up_cu;
+
+  gkyl_cu_memcpy(up_cu, &up_ho, sizeof(struct gkyl_proj_on_basis), GKYL_CU_MEMCPY_H2D);
+
+  // When the user does not provide a c2p mapping use the identity mapping,
+  // which must be assigned on the device to obtain a device function pointer.
+  if (up->c2p == 0)
+    proj_on_basis_set_c2p_identity_cu_ker<<<1,1>>>(up_cu);
+
+  return up_cu;
+}
+
+__global__ static void
+gkyl_proj_on_basis_advance_cu_ker(const struct gkyl_proj_on_basis *up, double tm,
+  struct gkyl_range update_range, struct gkyl_array* GKYL_RESTRICT fun_at_ords,
+  struct gkyl_array* GKYL_RESTRICT arr)
+{
+  int idx[GKYL_MAX_DIM];
+  double xc[GKYL_MAX_DIM], xmu[GKYL_MAX_DIM];
+
+  int num_basis = up->num_basis;
+  int tot_quad = up->tot_quad;
+  int num_ret_vals = up->num_ret_vals;
+  int ndim = up->grid.ndim;
+
+  const double* GKYL_RESTRICT weights = (const double*) up->weights->data;
+  const double* GKYL_RESTRICT basis_at_ords = (const double*) up->basis_at_ords->data;
+
+  for (unsigned long linc1 = threadIdx.x + blockIdx.x*blockDim.x;
+      linc1 < update_range.volume;
+      linc1 += gridDim.x*blockDim.x)
+  {
+    gkyl_sub_range_inv_idx(&update_range, linc1, idx);
+    gkyl_rect_grid_cell_center(&up->grid, idx, xc);
+
+    // Scratch space for the function evaluated at one quadrature node,
+    // one slot per cell (i.e. per linc1).
+    double *fq = (double*) gkyl_array_fetch(fun_at_ords, linc1);
+
+    long lidx = gkyl_range_idx(&update_range, idx);
+    double *f = (double*) gkyl_array_fetch(arr, lidx);
+
+    // Arrangement of f is as:
+    // c0[0], c0[1], ... c1[0], c1[1], ....
+    // where c0, c1, ... are components of f (num_ret_vals).
+    for (int k=0; k<num_ret_vals*num_basis; ++k) f[k] = 0.0;
+
+    for (int imu=0; imu<tot_quad; ++imu) {
+      proj_on_basis_log_to_comp(ndim, (const double*) gkyl_array_cfetch(up->ordinates, imu),
+        up->grid.dx, xc, xmu);
+      up->c2p(xmu, xmu, up->c2p_ctx);
+      up->eval(tm, xmu, fq, up->ctx);
+
+      long offset = 0;
+      for (int n=0; n<num_ret_vals; ++n) {
+        double tmp = weights[imu]*fq[n];
+        for (int k=0; k<num_basis; ++k)
+          f[offset+k] += tmp*basis_at_ords[k+num_basis*imu];
+        offset += num_basis;
+      }
+    }
+  }
+}
+
+void
+gkyl_proj_on_basis_advance_cu(const struct gkyl_proj_on_basis *up,
+  double tm, const struct gkyl_range *update_range, struct gkyl_array *arr)
+{
+  assert(gkyl_array_is_cu_dev(arr));
+
+  // Scratch storage for the function evaluated at a quadrature node, one slot
+  // per cell (allocated per call, like the CPU path allocates fun_at_ords).
+  struct gkyl_array *fun_at_ords = gkyl_array_cu_dev_new(GKYL_DOUBLE,
+    up->num_ret_vals, update_range->volume);
+
+  gkyl_proj_on_basis_advance_cu_ker<<<update_range->nblocks, update_range->nthreads>>>(
+    up->on_dev, tm, *update_range, fun_at_ords->on_dev, arr->on_dev);
+
+  gkyl_array_release(fun_at_ords);
+}
diff --git a/gyrokinetic/apps/gk_species_projection.c b/gyrokinetic/apps/gk_species_projection.c
index 843e0601e..fcc6763dd 100644
--- a/gyrokinetic/apps/gk_species_projection.c
+++ b/gyrokinetic/apps/gk_species_projection.c
@@ -67,23 +67,36 @@ func_gaussian(double t, const double* xn, double* GKYL_RESTRICT fout, void *ctx)
 }
 
 static void
-gk_species_projection_calc_proj_func(gkyl_gyrokinetic_app *app, struct gk_species *s, 
-  struct gk_proj *proj, struct gkyl_array *f, double tm) 
+gk_species_projection_multiply_jacobians(gkyl_gyrokinetic_app *app, struct gk_species *s,
+  struct gkyl_array *f)
 {
-  if (app->use_gpu) {
-    gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, proj->proj_host);
-    gkyl_array_copy(f, proj->proj_host);
-  }
-  else {
-    gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, f);
-  }
   // Multiply by the gyrocenter coord jacobian (bmag).
-  gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, f, 
-      app->gk_geom->geo_corn.bmag, f, &app->local, &s->local); 
-  // Multiply by the velocity-space jacobian. 
+  gkyl_dg_mul_conf_phase_op_range(&app->basis, &s->basis, f,
+      app->gk_geom->geo_corn.bmag, f, &app->local, &s->local);
+  // Multiply by the velocity-space jacobian.
   gkyl_array_scale_by_cell(f, s->vel_map->jacobvel);
 }
 
+static void
+gk_species_projection_calc_proj_func(gkyl_gyrokinetic_app *app, struct gk_species *s,
+  struct gk_proj *proj, struct gkyl_array *f, double tm)
+{
+  // The projection runs on the host or on the device (directly into f)
+  // depending on how proj_func was created.
+  gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, f);
+  gk_species_projection_multiply_jacobians(app, s, f);
+}
+
+static void
+gk_species_projection_calc_proj_func_host(gkyl_gyrokinetic_app *app, struct gk_species *s,
+  struct gk_proj *proj, struct gkyl_array *f, double tm)
+{
+  // Project on the host and copy to the device.
+  gkyl_proj_on_basis_advance(proj->proj_func, tm, &s->local, proj->proj_host);
+  gkyl_array_copy(f, proj->proj_host);
+  gk_species_projection_multiply_jacobians(app, s, f);
+}
+
 static void
 project_moment_if_needed(bool from_file, struct gkyl_proj_on_basis *proj_op,
   double tm, const struct gkyl_range *conf_range, struct gkyl_array *arr, double scale_fac)
@@ -430,22 +443,59 @@ gk_species_projection_init(struct gkyl_gyrokinetic_app *app, struct gk_species *
   proj->proj_on_basis_c2p_ctx.vel_map = s->vel_map;
   proj->proj_on_basis_c2p_ctx.pos_map = app->position_map;
   if (proj->proj_id == GKYL_PROJ_FUNC) {
-    proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
-        .grid = &s->grid,
-        .basis = &s->basis,
-        .qtype = GKYL_GAUSS_QUAD,
-        .num_quad = app->basis.poly_order+1,
-        .num_ret_vals = 1,
-        .eval = inp.func,
-        .ctx = inp.ctx_func,
-        .c2p_func = proj_on_basis_c2p_phase_func,
-        .c2p_func_ctx = &proj->proj_on_basis_c2p_ctx,
-      }
-    );
-    if (app->use_gpu)
-      proj->proj_host = mkarr(false, s->basis.num_basis, s->local_ext.volume);
+    proj->proj_host = 0;
+    proj->proj_on_basis_c2p_ctx_dev = 0;
+
+    bool proj_on_gpu = false;
+#ifdef GKYL_HAVE_CUDA
+    proj_on_gpu = app->use_gpu && (inp.func_on_dev != 0);
+    if (proj_on_gpu) {
+      // Project on the device with the user's device function. The c2p
+      // context must live on the device and hold the maps' device objects.
+      struct gk_proj_on_basis_c2p_func_ctx c2p_ctx_dev_ho = {
+        .cdim = app->cdim,
+        .vdim = s->local_vel.ndim,
+        .vel_map = s->vel_map->on_dev,
+        .pos_map = app->position_map->on_dev,
+      };
+      proj->proj_on_basis_c2p_ctx_dev = gkyl_cu_malloc(sizeof(struct gk_proj_on_basis_c2p_func_ctx));
+      gkyl_cu_memcpy(proj->proj_on_basis_c2p_ctx_dev, &c2p_ctx_dev_ho,
+        sizeof(struct gk_proj_on_basis_c2p_func_ctx), GKYL_CU_MEMCPY_H2D);
+
+      proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+          .grid = &s->grid,
+          .basis = &s->basis,
+          .qtype = GKYL_GAUSS_QUAD,
+          .num_quad = app->basis.poly_order+1,
+          .num_ret_vals = 1,
+          .eval = inp.func_on_dev,
+          .ctx = inp.ctx_func_on_dev,
+          .c2p_func = gk_species_projection_c2p_phase_func_cu_dev_ptr(),
+          .c2p_func_ctx = proj->proj_on_basis_c2p_ctx_dev,
+          .use_gpu = true,
+        }
+      );
+    }
+#endif
+    if (!proj_on_gpu) {
+      proj->proj_func = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
+          .grid = &s->grid,
+          .basis = &s->basis,
+          .qtype = GKYL_GAUSS_QUAD,
+          .num_quad = app->basis.poly_order+1,
+          .num_ret_vals = 1,
+          .eval = inp.func,
+          .ctx = inp.ctx_func,
+          .c2p_func = proj_on_basis_c2p_phase_func,
+          .c2p_func_ctx = &proj->proj_on_basis_c2p_ctx,
+        }
+      );
+      if (app->use_gpu)
+        proj->proj_host = mkarr(false, s->basis.num_basis, s->local_ext.volume);
+    }
 
-    proj->projection_calc = gk_species_projection_calc_proj_func;
+    proj->projection_calc = (app->use_gpu && !proj_on_gpu) ?
+      gk_species_projection_calc_proj_func_host : gk_species_projection_calc_proj_func;
     proj->moms_correct = gk_species_projection_correct_all_moms_none;
   }
   else {
@@ -478,9 +528,14 @@ gk_species_projection_release(const struct gkyl_gyrokinetic_app *app, const stru
 {
   if (proj->proj_id == GKYL_PROJ_FUNC) {
     gkyl_proj_on_basis_release(proj->proj_func);
-    if (app->use_gpu) {
+    if (proj->proj_host) {
       gkyl_array_release(proj->proj_host);
     }
+#ifdef GKYL_HAVE_CUDA
+    if (proj->proj_on_basis_c2p_ctx_dev) {
+      gkyl_cu_free(proj->proj_on_basis_c2p_ctx_dev);
+    }
+#endif
   }
   else if (proj->proj_id == GKYL_PROJ_MAXWELLIAN_PRIM || proj->proj_id == GKYL_PROJ_BIMAXWELLIAN) { 
     gkyl_array_release(proj->dens);
diff --git a/gyrokinetic/apps/gk_species_projection_cu.cu b/gyrokinetic/apps/gk_species_projection_cu.cu
new file mode 100644
index 000000000..a493517a4
--- /dev/null
+++ b/gyrokinetic/apps/gk_species_projection_cu.cu
@@ -0,0 +1,40 @@
+/* -*- c++ -*- */
+
+extern "C" {
+#include <gkyl_gk_proj_on_basis_c2p_priv.h>
+#include <gkyl_position_map.h>
+#include <gkyl_proj_on_basis.h>
+#include <gkyl_util.h>
+#include <gkyl_velocity_map.h>
+}
+
+// Comp. to phys. mapping for phase-space projections (position map in
+// configuration space, velocity map in velocity space), callable on the
+// device. The context must be a GPU-resident gk_proj_on_basis_c2p_func_ctx
+// whose pos_map/vel_map members are the maps' device (on_dev) objects.
+GKYL_CU_D static void
+gk_proj_on_basis_c2p_phase_func_cu(const double *xcomp, double *xphys, void *ctx)
+{
+  struct gk_proj_on_basis_c2p_func_ctx *c2p_ctx = (struct gk_proj_on_basis_c2p_func_ctx *) ctx;
+  int cdim = c2p_ctx->cdim; // Assumes update range is a phase range.
+  gkyl_position_map_eval_mc2nu(c2p_ctx->pos_map, xcomp, xphys);
+  gkyl_velocity_map_eval_c2p_dev(c2p_ctx->vel_map, &xcomp[cdim], &xphys[cdim]);
+}
+
+__global__ static void
+gk_proj_on_basis_c2p_phase_func_set_cu_ker(proj_on_basis_c2p_t *fptr_d)
+{
+  // Assigned in device code so the function pointer is a device address.
+  *fptr_d = gk_proj_on_basis_c2p_phase_func_cu;
+}
+
+extern "C" proj_on_basis_c2p_t
+gk_species_projection_c2p_phase_func_cu_dev_ptr(void)
+{
+  proj_on_basis_c2p_t *fptr_d, fptr_ho;
+  checkCuda(cudaMalloc(&fptr_d, sizeof(proj_on_basis_c2p_t)));
+  gk_proj_on_basis_c2p_phase_func_set_cu_ker<<<1,1>>>(fptr_d);
+  checkCuda(cudaMemcpy(&fptr_ho, fptr_d, sizeof(proj_on_basis_c2p_t), cudaMemcpyDeviceToHost));
+  checkCuda(cudaFree(fptr_d));
+  return fptr_ho;
+}
diff --git a/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h b/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h
new file mode 100644
index 000000000..c2faa069e
--- /dev/null
+++ b/gyrokinetic/apps/gkyl_gk_proj_on_basis_c2p_priv.h
@@ -0,0 +1,29 @@
+// Private header for use in the Gyrokinetic app: do not include in
+// user-facing header files! Context (and device helpers) for the
+// computational-to-physical coordinate mappings passed to proj_on_basis.
+#pragma once
+
+#include <gkyl_position_map.h>
+#include <gkyl_proj_on_basis.h>
+#include <gkyl_velocity_map.h>
+
+// Context for c2p function passed to proj_on_basis.
+struct gk_proj_on_basis_c2p_func_ctx {
+  int cdim, vdim;
+  struct gkyl_position_map *pos_map;
+  struct gkyl_velocity_map *vel_map;
+};
+
+#ifdef GKYL_HAVE_CUDA
+
+/**
+ * Return the device address of the phase-space c2p function (position map in
+ * configuration space, velocity map in velocity space), for projecting with
+ * proj_on_basis on the GPU. The context passed along with it must be a
+ * GPU-resident gk_proj_on_basis_c2p_func_ctx whose pos_map/vel_map members
+ * are the maps' device (on_dev) objects. Defined in
+ * gk_species_projection_cu.cu.
+ */
+proj_on_basis_c2p_t gk_species_projection_c2p_phase_func_cu_dev_ptr(void);
+
+#endif
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 1d137b096..1c44cd349 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -38,8 +38,18 @@ struct gkyl_gyrokinetic_projection {
       // Distribution function to project and its context. For fluid neutrals
       // this function returns mass density, momentum density, and total
       // energy density.
-      void (*func)(double t, const double *xn, double *fout, void *ctx); 
-      void *ctx_func; 
+      void (*func)(double t, const double *xn, double *fout, void *ctx);
+      void *ctx_func;
+
+      // Optionally, when running on GPUs, the device address of a
+      // device-callable version of 'func' (defined in a CUDA-compiled
+      // translation unit and obtained with GKYL_DEFINE_CU_DEV_FUNC_GETTER)
+      // and a device-resident context. When provided, the projection is
+      // performed on the GPU; the user must place the context on the GPU
+      // themselves as the app has no way to perform that copy. If absent,
+      // the function is projected on the host and copied to the GPU.
+      void (*func_on_dev)(double t, const double *xn, double *fout, void *ctx);
+      void *ctx_func_on_dev;
     };
     struct {
       // For Maxwellians (or BiMaxwellians), specify density, parallel speed
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index e8d095590..469a1abd6 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -724,12 +724,9 @@ struct gk_scaling {
     struct gk_scaling *sca, int ridx, double tm, int frame);
 };
 
-// Context for c2p function passed to proj_on_basis.
-struct gk_proj_on_basis_c2p_func_ctx {
-  int cdim, vdim;
-  struct gkyl_position_map *pos_map;
-  struct gkyl_velocity_map *vel_map;
-};
+// Context for c2p function passed to proj_on_basis is defined in
+// gkyl_gk_proj_on_basis_c2p_priv.h (shared with gk_species_projection_cu.cu).
+#include <gkyl_gk_proj_on_basis_c2p_priv.h>
 
 struct gk_proj {
   enum gkyl_projection_id proj_id; // Type of projection.
@@ -739,6 +736,8 @@ struct gk_proj {
     struct {
       struct gkyl_proj_on_basis *proj_func; // Projection operator for specified function.
       struct gkyl_array *proj_host; // Array for projection on host-side if running on GPUs.
+      struct gk_proj_on_basis_c2p_func_ctx *proj_on_basis_c2p_ctx_dev; // Device copy of the
+                                                                       // c2p context (GPU projection).
     };
     // Maxwellian and Bi-Maxwellian projection from primitive moments.
     struct {
diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index 5cf1ff2c9..7b2b5e493 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -417,6 +417,14 @@ gkyl_gyrokinetic_app_new_geom(struct gkyl_gk *gk)
 
   gkyl_position_map_set_mc2nu(app->position_map, app->gk_geom->geo_corn.mc2nu_pos);
 
+#ifdef GKYL_HAVE_CUDA
+  if (app->use_gpu) {
+    // Device copy of the position map, used e.g. by projections that
+    // evaluate the c2p mapping inside device kernels.
+    gkyl_position_map_make_cu_dev(app->position_map);
+  }
+#endif
+
   const struct gkyl_dg_geom_inp dg_geom_inp = {
     .grid = &app->grid,
     .range = &app->local_ext,
diff --git a/gyrokinetic/zero/gkyl_position_map.h b/gyrokinetic/zero/gkyl_position_map.h
index 8d88742e8..e3493566d 100644
--- a/gyrokinetic/zero/gkyl_position_map.h
+++ b/gyrokinetic/zero/gkyl_position_map.h
@@ -6,6 +6,9 @@
 #include <gkyl_array.h>
 #include <gkyl_comm.h>
 #include <gkyl_comm_io.h>
+#include <gkyl_util.h>
+
+#include <math.h>
 
 enum gkyl_position_map_id {
   GKYL_PMAP_USER_INPUT = 0, // Function projection. User specified. Default
@@ -62,7 +65,11 @@ struct gkyl_position_map {
   // Stuff for constant B mapping
   struct gkyl_bmag_ctx *bmag_ctx; // Context for magnetic field calculation
   struct gkyl_position_map_const_B_ctx *constB_ctx; // Context for constant B mapping
-  struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping 
+  struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping
+
+  struct gkyl_array *mc2nu_dev; // Device mirror of mc2nu (only with a device copy, see
+                                // gkyl_position_map_make_cu_dev; kept in sync by set_mc2nu).
+  struct gkyl_position_map *on_dev; // Device copy of itself (points to itself on CPU).
 };
 
 struct gkyl_position_map_const_B_ctx {
@@ -175,16 +182,54 @@ void
 gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, 
     double zcenter, double w, double psisep);
 
+/**
+ * Create a device (GPU) copy of the position map object, stored in
+ * gpm->on_dev, for use inside device kernels with
+ * gkyl_position_map_eval_mc2nu. Call it after the geometry has set the
+ * mc2nu array (gkyl_position_map_set_mc2nu keeps the device mirror in sync
+ * afterwards). Only available in GPU builds.
+ *
+ * @param gpm Position map object.
+ */
+void
+gkyl_position_map_make_cu_dev(struct gkyl_position_map* gpm);
+
 /**
  * Evaluate the position mapping at a specific computational (position) coordinate.
- * NOTE: done on the host.
+ * Callable on the host with the host object, or inside device kernels with the
+ * device object (gpm->on_dev). On GPU builds do not call it on the host with
+ * the device object or vice-versa.
  *
  * @param gpm Gkyl position map object.
- * @param xc Computational position coordinates.
- * @param xnu Resulting non-uniform position coordinates.
+ * @param x_comp Computational position coordinates.
+ * @param x_fa Resulting non-uniform position coordinates.
  */
-void
-gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *xc, double *xnu);
+GKYL_CU_DH static inline void
+gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *x_comp, double *x_fa)
+{
+  int cidx[GKYL_MAX_CDIM];
+  for (int i=0; i<gpm->grid.ndim; i++) {
+    int idxtemp = gpm->global.lower[i] + (int) floor((x_comp[i] - (gpm->grid.lower[i]) )/gpm->grid.dx[i]);
+    idxtemp = GKYL_MAX2(GKYL_MIN2(idxtemp, gpm->local.upper[i]), gpm->local.lower[i]);
+    cidx[i] = idxtemp;
+  }
+  long lidx = gkyl_range_idx(&gpm->local, cidx);
+  const double *pmap_coeffs = (const double *) gkyl_array_cfetch(gpm->mc2nu, lidx);
+  double cxc[GKYL_MAX_CDIM];
+  double x_log[GKYL_MAX_CDIM];
+  gkyl_rect_grid_cell_center(&gpm->grid, cidx, cxc);
+  for (int i=0; i<gpm->grid.ndim; i++) {
+    x_log[i] = (x_comp[i]-cxc[i])/(gpm->grid.dx[i]*0.5);
+  }
+  double xyz_fa[3];
+  for (int i=0; i<3; i++) {
+    xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i*gpm->basis.num_basis]);
+  }
+  for (int i=0; i<gpm->grid.ndim; i++) {
+    x_fa[i] = xyz_fa[i];
+  }
+  x_fa[gpm->grid.ndim-1] = xyz_fa[2];
+}
 
 /**
  * Evaluate the slope of the position mapping at a specific computational (position) coordinate.
diff --git a/gyrokinetic/zero/position_map.c b/gyrokinetic/zero/position_map.c
index 47b0e9523..9a2e70a6c 100644
--- a/gyrokinetic/zero/position_map.c
+++ b/gyrokinetic/zero/position_map.c
@@ -9,6 +9,7 @@
 
 #include <float.h>
 #include <assert.h>
+#include <stddef.h>
 
 // Remove with the print statements at the bottom
 #include <gkyl_util.h>
@@ -39,6 +40,8 @@ gkyl_position_map_null_new()
   gpm->xpt_ctx = gkyl_malloc(sizeof(struct gkyl_position_map_xpt_ctx));
   gpm->bmag_ctx = gkyl_malloc(sizeof(struct gkyl_bmag_ctx));
   gpm->bmag_ctx->bmag = gkyl_array_new(GKYL_DOUBLE, 1, 1);
+  gpm->mc2nu_dev = 0;
+  gpm->on_dev = gpm; // On the CPU the object points to itself.
   gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free);
   
   for (int i = 0; i < 3; i++){
@@ -161,6 +164,8 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g
   gpm->basis = basis;
   gpm->cdim = grid.ndim; 
   gpm->mc2nu = gkyl_array_new(GKYL_DOUBLE, 3*gpm->basis.num_basis, gpm->local_ext.volume);
+  gpm->mc2nu_dev = 0;
+  gpm->on_dev = gpm; // On the CPU the object points to itself.
   gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free);
 
   struct gkyl_position_map *gpm_out = gpm;
@@ -171,6 +176,8 @@ void
 gkyl_position_map_set_mc2nu(struct gkyl_position_map* gpm, struct gkyl_array* mc2nu)
 {
   gkyl_array_copy(gpm->mc2nu, mc2nu);
+  if (gpm->mc2nu_dev)
+    gkyl_array_copy(gpm->mc2nu_dev, gpm->mc2nu);
 }
 
 void
@@ -225,31 +232,53 @@ gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, do
   }
 }
 
-void 
-gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *x_comp, double *x_fa)
+void
+gkyl_position_map_make_cu_dev(struct gkyl_position_map* gpm)
 {
-  int cidx[GKYL_MAX_CDIM];
-  for(int i = 0; i < gpm->grid.ndim; i++){
-    int idxtemp = gpm->global.lower[i] + (int) floor((x_comp[i] - (gpm->grid.lower[i]) )/gpm->grid.dx[i]);
-    idxtemp = GKYL_MAX2(GKYL_MIN2(idxtemp, gpm->local.upper[i]), gpm->local.lower[i]);
-    cidx[i] = idxtemp;
-  }
-  long lidx = gkyl_range_idx(&gpm->local, cidx);
-  const double *pmap_coeffs = gkyl_array_cfetch(gpm->mc2nu, lidx);
-  double cxc[gpm->grid.ndim];
-  double x_log[gpm->grid.ndim];
-  gkyl_rect_grid_cell_center(&gpm->grid, cidx, cxc);
-  for(int i = 0; i < gpm->grid.ndim; i++){
-    x_log[i] = (x_comp[i]-cxc[i])/(gpm->grid.dx[i]*0.5);
-  }
-  double xyz_fa[3];
-  for(int i = 0; i < 3; i++){
-    xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i*gpm->basis.num_basis]);
+#ifdef GKYL_HAVE_CUDA
+  // Device mirror of the mapping's DG coefficients, kept in sync by set_mc2nu.
+  gpm->mc2nu_dev = gkyl_array_cu_dev_new(GKYL_DOUBLE, gpm->mc2nu->ncomp, gpm->mc2nu->size);
+  gkyl_array_copy(gpm->mc2nu_dev, gpm->mc2nu);
+
+  // Clone with device pointers; host-only members are nulled.
+  struct gkyl_position_map gpm_ho = *gpm;
+  gpm_ho.mc2nu = gpm->mc2nu_dev->on_dev;
+  for (int i=0; i<3; i++) {
+    gpm_ho.maps[i] = 0;
+    gpm_ho.map_derivs[i] = 0;
+    gpm_ho.ctxs[i] = 0;
   }
-  for (int i=0; i<gpm->grid.ndim; i++) {
-    x_fa[i] = xyz_fa[i];
+  gpm_ho.bmag_ctx = 0;
+  gpm_ho.constB_ctx = 0;
+  gpm_ho.xpt_ctx = 0;
+
+  struct gkyl_position_map *gpm_cu = gkyl_cu_malloc(sizeof(struct gkyl_position_map));
+  gpm_ho.on_dev = gpm_cu;
+  gkyl_cu_memcpy(gpm_cu, &gpm_ho, sizeof(struct gkyl_position_map), GKYL_CU_MEMCPY_H2D);
+
+  // Overwrite the basis in the clone with one whose function pointers are
+  // device addresses (the temporary device container can be released after
+  // copying its contents; the code addresses it holds remain valid).
+  struct gkyl_basis *basis_cu;
+  switch (gpm->basis.b_type) {
+    case GKYL_BASIS_MODAL_SERENDIPITY:
+      basis_cu = gkyl_cart_modal_serendip_cu_dev_new(gpm->basis.ndim, gpm->basis.poly_order);
+      break;
+    case GKYL_BASIS_MODAL_TENSOR:
+      basis_cu = gkyl_cart_modal_tensor_cu_dev_new(gpm->basis.ndim, gpm->basis.poly_order);
+      break;
+    default:
+      assert(false);
+      break;
   }
-  x_fa[gpm->grid.ndim-1] = xyz_fa[2];
+  gkyl_cu_memcpy((char *)gpm_cu + offsetof(struct gkyl_position_map, basis), basis_cu,
+    sizeof(struct gkyl_basis), GKYL_CU_MEMCPY_D2D);
+  gkyl_cart_modal_basis_release_cu(basis_cu);
+
+  gpm->on_dev = gpm_cu;
+#else
+  assert(false);
+#endif
 }
 
 void
@@ -364,6 +393,10 @@ gkyl_position_map_free(const struct gkyl_ref_count *ref)
   struct gkyl_position_map *gpm = container_of(ref, struct gkyl_position_map, ref_count);
   gkyl_array_release(gpm->mc2nu);
   gkyl_array_release(gpm->bmag_ctx->bmag);
+  if (gpm->mc2nu_dev)
+    gkyl_array_release(gpm->mc2nu_dev);
+  if (gpm->on_dev != gpm)
+    gkyl_cu_free(gpm->on_dev);
   if (gpm->to_optimize == true)
   {
     gkyl_free(gpm->constB_ctx->theta_extrema);
diff --git a/vlasov/zero/gkyl_velocity_map.h b/vlasov/zero/gkyl_velocity_map.h
index b7054bfa4..17ce622f9 100644
--- a/vlasov/zero/gkyl_velocity_map.h
+++ b/vlasov/zero/gkyl_velocity_map.h
@@ -6,6 +6,9 @@
 #include <gkyl_array.h>
 #include <gkyl_comm.h>
 #include <gkyl_comm_io.h>
+#include <gkyl_util.h>
+
+#include <math.h>
 
 typedef void (*mapc2p_t)(double t, const double *zc, double *vp, void *ctx);
 
@@ -115,6 +118,47 @@ gkyl_velocity_map_reduce_dv_range(const struct gkyl_velocity_map* gvm, enum gkyl
 void
 gkyl_velocity_map_eval_c2p(const struct gkyl_velocity_map* gvm, const double *zc, double *vp);
 
+/**
+ * Evaluate the velocity mapping at a specific computational (velocity)
+ * coordinate. Unlike gkyl_velocity_map_eval_c2p, this variant uses the vmap
+ * and vmap_basis members (not their _ho host copies), so it is callable
+ * inside device kernels when given the device object (gvm->on_dev). On GPU
+ * builds do not call it on the host with the host object, whose vmap data
+ * lives on the device.
+ *
+ * @param gvm Velocity map object (on_dev object inside kernels).
+ * @param zc Computational velocity coordinates.
+ * @param vp Resulting physical velocity coordinates.
+ */
+GKYL_CU_DH static inline void
+gkyl_velocity_map_eval_c2p_dev(const struct gkyl_velocity_map* gvm, const double *zc, double *vp)
+{
+  // Find the index of the cell containing zc.
+  int idx_zc[GKYL_MAX_VDIM];
+  for (int d=0; d<gvm->local_ext_vel.ndim; d++) {
+    int idx = gvm->local_ext_vel.lower[d] + (int) floor((zc[d] - (gvm->grid_vel.lower[d]) )/gvm->grid_vel.dx[d]);
+    // Bound idx to the range in the grid. If it falls outside of that is due
+    // to floating point arithmetic, or due to an error in the code.
+    idx = GKYL_MIN2(idx, gvm->local_ext_vel.upper[d]);
+    idx = GKYL_MAX2(idx, gvm->local_ext_vel.lower[d]);
+    idx_zc[d] = idx;
+  }
+
+  // Fetch DG coefficients of the velocity map in idx_zc.
+  long lidx_zc = gkyl_range_idx(&gvm->local_ext_vel, idx_zc);
+  const double *vmap_c = (const double *) gkyl_array_cfetch(gvm->vmap, lidx_zc);
+
+  double zc_cc[GKYL_MAX_VDIM];
+  gkyl_rect_grid_cell_center(&gvm->grid_vel, idx_zc, zc_cc);
+
+  for (int d=0; d<gvm->local_ext_vel.ndim; d++) {
+    // Convert computational to logical coord.
+    double zlog[] = {(zc[d] - zc_cc[d]) / (0.5*gvm->grid_vel.dx[d])};
+    // Evaluate vmap expansion at logical coord.
+    vp[d] = gvm->vmap_basis->eval_expand(zlog, &vmap_c[d*gvm->vmap_basis->num_basis]);
+  }
+}
+
 /**
  * Indicate if this velocity map object is allocated on the GPU.
  *