array_dg_find_peaks.c

#include <assert.h>
#include <float.h>
#include <string.h>

#include <gkyl_alloc.h>
#include <gkyl_alloc_flags_priv.h>
#include <gkyl_array.h>
#include <gkyl_array_dg_find_peaks.h>
#include <gkyl_array_dg_find_peaks_priv.h>
#include <gkyl_nodal_ops.h>

/**
 * Scan along the search direction at a fixed preserved-direction coordinate
 * to count the number of peaks and determine their types.
 */
static void
count_peaks_along_dir(const struct gkyl_array_dg_find_peaks *up, const struct gkyl_array *in_ho,
  int preserved_idx, int *num_peaks_out, enum gkyl_peak_type *peak_types_out)
{
  int ndim = up->grid.ndim;
  int search_dir = up->search_dir;

  int total_nodes_search = up->total_nodes_search;

  // Use pre-allocated search buffers from the struct.
  double *vals = up->search_vals;
  double *coords = up->search_coords;
  for (int i = 0; i < total_nodes_search; i++) {
    vals[i] = 0.0;
    coords[i] = 0.0;
  }

  // Iterate along cells in search direction and collect nodal values.
  for (int cell_idx = up->range.lower[search_dir];
    cell_idx <= up->range.upper[search_dir];
    cell_idx++) {
    // Build index array for this cell.
    int idx[GKYL_MAX_DIM];
    if (ndim == 1) {
      idx[0] = cell_idx;
    }
    else {
      int preserved_dir = (search_dir == 0) ? 1 : 0;
      idx[preserved_dir] = preserved_idx;
      idx[search_dir] = cell_idx;
    }

    long linidx = gkyl_range_idx(&up->range, idx);
    const double *f_d = gkyl_array_cfetch(in_ho, linidx);

    double xc[GKYL_MAX_DIM];
    gkyl_rect_grid_cell_center(&up->grid, idx, xc);

    // Evaluate at each node in this cell.
    for (int n = 0; n < up->basis.num_basis; n++) {
      const double *nod_log = gkyl_array_cfetch(up->nodes, n);

      // Determine node offset in search direction.
      int node_offset = (nod_log[search_dir] < 0) ? 0 : 1;

      int cell_local = cell_idx - up->range.lower[search_dir];

      int search_node_idx = cell_local + node_offset;

      double val = up->basis.eval_expand(nod_log, f_d);
      double nod_phys[GKYL_MAX_DIM];
      dg_find_peaks_log_to_comp(ndim, nod_log, up->grid.dx, xc, nod_phys);

      // Only store if this is the first time we see this search node
      // (avoid duplicates at cell boundaries).
      if (vals[search_node_idx] == 0.0 && coords[search_node_idx] == 0.0) {
        vals[search_node_idx] = val;
        coords[search_node_idx] = nod_phys[search_dir];
      }
    }
  }

  // Now scan the values to find peaks.
  // A peak is: EDGE_LO at index 0, EDGE_HI at last index, LOCAL_MAX/MIN in between.
  int num_peaks = 0;

  // Always add lower edge.
  peak_types_out[num_peaks++] = GKYL_PEAK_EDGE_LO;

  // Scan for local maxima and minima (indices 1 to total_nodes_search-2).
  for (int i = 1; i < total_nodes_search - 1; i++) {
    double prev = vals[i - 1];
    double curr = vals[i];
    double next = vals[i + 1];

    if (curr > prev && curr > next) {
      // Local maximum.
      assert(num_peaks < GKYL_DG_FIND_PEAKS_MAX);
      peak_types_out[num_peaks++] = GKYL_PEAK_LOCAL_MAX;
    }
    else if (curr < prev && curr < next) {
      // Local minimum.
      assert(num_peaks < GKYL_DG_FIND_PEAKS_MAX);
      peak_types_out[num_peaks++] = GKYL_PEAK_LOCAL_MIN;
    }
  }

  // Always add upper edge.
  assert(num_peaks < GKYL_DG_FIND_PEAKS_MAX);
  peak_types_out[num_peaks++] = GKYL_PEAK_EDGE_HI;

  *num_peaks_out = num_peaks;
}

/**
 * Find all peaks along the search direction for a given preserved-direction
 * node index, storing results in the nodal arrays.
 */
static void
find_peaks_for_preserved_node(struct gkyl_array_dg_find_peaks *up, const struct gkyl_array *in_ho,
  int preserved_node_idx)
{
  int ndim = up->grid.ndim;
  int search_dir = up->search_dir;

  int total_nodes_search = up->total_nodes_search;

  // Use pre-allocated search buffers from the struct.
  double *vals = up->search_vals;
  double *coords = up->search_coords;
  bool *visited = up->search_visited;
  for (int i = 0; i < total_nodes_search; i++) {
    vals[i] = 0.0;
    coords[i] = 0.0;
    visited[i] = false;
  }

  // For 2D, determine the preserved direction cell index from the node index.
  int preserved_dir = (ndim == 1) ? -1 : ((search_dir == 0) ? 1 : 0);

  // Iterate along cells in search direction and collect nodal values.
  for (int cell_idx = up->range.lower[search_dir];
    cell_idx <= up->range.upper[search_dir];
    cell_idx++) {
    // For 2D, we need to iterate over cells in the preserved direction that
    // contribute to this preserved node index.
    int pres_cell_start, pres_cell_end;
    if (ndim == 1) {
      pres_cell_start = 0;
      pres_cell_end = 0;
    }
    else {
      // Determine which cells contribute to this preserved node.
      // Node i is shared by cells i and i+1 (0-indexed from lower).
      // preserved_node_idx 0 is only in cell lower[preserved_dir].
      // preserved_node_idx N is only in cell upper[preserved_dir].
      if (preserved_node_idx == 0) {
        pres_cell_start = up->range.lower[preserved_dir];
        pres_cell_end = up->range.lower[preserved_dir];
      }
      else if (preserved_node_idx == up->out_nrange.upper[0]) {
        pres_cell_start = up->range.upper[preserved_dir];
        pres_cell_end = up->range.upper[preserved_dir];
      }
      else {
        pres_cell_start = up->range.lower[preserved_dir] + preserved_node_idx - 1;
        pres_cell_end = pres_cell_start + 1;
        if (pres_cell_end > up->range.upper[preserved_dir])
          pres_cell_end = up->range.upper[preserved_dir];
      }
    }

    for (int pres_cell = pres_cell_start; pres_cell <= pres_cell_end; pres_cell++) {
      // Build index array for this cell.
      int idx[GKYL_MAX_DIM];
      if (ndim == 1) {
        idx[0] = cell_idx;
      }
      else {
        idx[preserved_dir] = pres_cell;
        idx[search_dir] = cell_idx;
      }

      long linidx = gkyl_range_idx(&up->range, idx);
      const double *f_d = gkyl_array_cfetch(in_ho, linidx);

      double xc[GKYL_MAX_DIM];
      gkyl_rect_grid_cell_center(&up->grid, idx, xc);

      // Evaluate at each node in this cell.
      for (int n = 0; n < up->basis.num_basis; n++) {
        const double *nod_log = gkyl_array_cfetch(up->nodes, n);

        // Check if this node corresponds to our preserved node index.
        if (ndim > 1) {
          int pres_node_offset = (nod_log[preserved_dir] < 0) ? 0 : 1;
          int pres_cell_local = pres_cell - up->range.lower[preserved_dir];
          int this_pres_node = pres_cell_local + pres_node_offset;
          if (this_pres_node != preserved_node_idx)
            continue;
        }

        // Determine node offset in search direction.
        int search_node_offset = (nod_log[search_dir] < 0) ? 0 : 1;

        int cell_local = cell_idx - up->range.lower[search_dir];
        int search_node_idx = cell_local + search_node_offset;

        if (!visited[search_node_idx]) {
          double val = up->basis.eval_expand(nod_log, f_d); // GPU error here
          double nod_phys[GKYL_MAX_DIM];
          dg_find_peaks_log_to_comp(ndim, nod_log, up->grid.dx, xc, nod_phys);

          vals[search_node_idx] = val;
          coords[search_node_idx] = nod_phys[search_dir];
          visited[search_node_idx] = true;
        }
      }
    }
  }

  // Now extract peaks based on peak_types.
  int peak_idx = 0;

  // EDGE_LO is always first peak at index 0.
  if (up->peak_types[peak_idx] == GKYL_PEAK_EDGE_LO) {
    double *val_n = gkyl_array_fetch(up->out_vals_nodal[peak_idx], preserved_node_idx);
    double *coord_n = gkyl_array_fetch(up->out_coords_nodal[peak_idx], preserved_node_idx);
    val_n[0] = vals[0];
    coord_n[0] = coords[0];
    peak_idx++;
  }

  // Find local maxima and minima.
  for (int i = 1; i < total_nodes_search - 1 && peak_idx < up->num_peaks - 1; i++) {
    double prev = vals[i - 1];
    double curr = vals[i];
    double next = vals[i + 1];

    bool is_max = (curr > prev && curr > next);
    bool is_min = (curr < prev && curr < next);

    if ((is_max && up->peak_types[peak_idx] == GKYL_PEAK_LOCAL_MAX) ||
      (is_min && up->peak_types[peak_idx] == GKYL_PEAK_LOCAL_MIN)) {
      double *val_n = gkyl_array_fetch(up->out_vals_nodal[peak_idx], preserved_node_idx);
      double *coord_n = gkyl_array_fetch(up->out_coords_nodal[peak_idx], preserved_node_idx);
      val_n[0] = curr;
      coord_n[0] = coords[i];
      peak_idx++;
    }
  }

  // EDGE_HI is always last peak.
  if (peak_idx < up->num_peaks && up->peak_types[peak_idx] == GKYL_PEAK_EDGE_HI) {
    double *val_n = gkyl_array_fetch(up->out_vals_nodal[peak_idx], preserved_node_idx);
    double *coord_n = gkyl_array_fetch(up->out_coords_nodal[peak_idx], preserved_node_idx);
    val_n[0] = vals[total_nodes_search - 1];
    coord_n[0] = coords[total_nodes_search - 1];
  }
}

/**
 * Evaluate an input array at peak locations for a given preserved-direction
 * node index, storing results in the nodal output arrays.
 */
static void
eval_array_at_peaks_for_preserved_node(struct gkyl_array_dg_find_peaks *up,
  const struct gkyl_array *in_ho, int preserved_node_idx, struct gkyl_array **out_vals_nodal,
  int peak_idx)
{
  int ndim = up->grid.ndim;
  int search_dir = up->search_dir;
  int preserved_dir = (ndim == 1) ? -1 : ((search_dir == 0) ? 1 : 0);

  // Get the peak coordinate that was found during find_peaks.
  const double *peak_coord_n = gkyl_array_cfetch(up->out_coords_nodal[peak_idx],
    preserved_node_idx);
  double peak_coord_search = peak_coord_n[0];

  // Find the cell containing this coordinate in the search direction.
  // We need to build a point coordinate to pass to find_cell.
  double point[GKYL_MAX_DIM];
  int known_idx[GKYL_MAX_DIM];
  int cell_idx[GKYL_MAX_DIM];

  for (int d = 0; d < ndim; d++) {
    if (d == search_dir) {
      point[d] = peak_coord_search;
      known_idx[d] = -1; // Not known
    }
    else {
      // Use dummy value - we'll specify known_idx.
      point[d] = 0.0;
      known_idx[d] = -1;
    }
  }

  // If 2D, we need to determine preserved direction cell from preserved_node_idx.
  // For p=1 with N cells (1-based indexing), nodal points map as:
  // Node 0 -> cell 1, logical coord -1 (left edge of first cell)
  // Node k (1 <= k <= N) -> cell k, logical coord +1 (right edge of cell k)
  // This ensures proper continuity at shared cell boundaries.
  if (ndim > 1) {
    int pres_cell;
    if (preserved_node_idx == 0) {
      // First node: evaluate at left edge of first cell.
      pres_cell = up->range.lower[preserved_dir];
    }
    else {
      // All other nodes (1 to N): evaluate at right edge of cell with index = node_idx.
      // Clamp to upper bound for safety.
      pres_cell = up->range.lower[preserved_dir] + preserved_node_idx - 1;
      if (pres_cell > up->range.upper[preserved_dir]) {
        pres_cell = up->range.upper[preserved_dir];
      }
    }
    known_idx[preserved_dir] = pres_cell;

    // Set the coordinate in preserved direction to the cell center.
    int pres_cell_idx[GKYL_MAX_DIM];
    for (int d = 0; d < ndim; d++) {
      pres_cell_idx[d] = (d == preserved_dir) ? pres_cell : 1;
    }
    double xc_pres[GKYL_MAX_DIM];
    gkyl_rect_grid_cell_center(&up->grid, pres_cell_idx, xc_pres);
    point[preserved_dir] = xc_pres[preserved_dir];
  }

  gkyl_rect_grid_find_cell(&up->grid, point, true, known_idx, cell_idx);

  // Clamp cell_idx to interior range (avoid ghost cells).
  for (int d = 0; d < up->grid.ndim; d++) {
    if (cell_idx[d] < up->range.lower[d]) {
      cell_idx[d] = up->range.lower[d];
    }
    if (cell_idx[d] > up->range.upper[d]) {
      cell_idx[d] = up->range.upper[d];
    }
  }

  // Get the DG coefficients at this cell.
  long linidx = gkyl_range_idx(&up->range, cell_idx);
  const double *f_d = gkyl_array_cfetch(in_ho, linidx);

  // Get cell center.
  double xc[GKYL_MAX_DIM];
  gkyl_rect_grid_cell_center(&up->grid, cell_idx, xc);

  // Convert peak coordinate to logical space.
  double nod_log[GKYL_MAX_DIM];
  for (int d = 0; d < ndim; d++) {
    if (d == search_dir) {
      // Convert physical coordinate to logical [-1, 1].
      nod_log[d] = 2.0 * (peak_coord_search - xc[d]) / up->grid.dx[d];
    }
    else if (ndim > 1) {
      // In preserved direction, use the node position in the cell.
      // Node 0 is at left edge (-1), all others at right edge (+1).
      nod_log[d] = (preserved_node_idx == 0) ? -1.0 : 1.0;
    }
  }

  // Evaluate the DG expansion at this logical coordinate.
  double val = up->basis.eval_expand(nod_log, f_d);

  // Store the result.
  double *val_n = gkyl_array_fetch(out_vals_nodal[peak_idx], preserved_node_idx);
  val_n[0] = val;
}

struct gkyl_array_dg_find_peaks*
gkyl_array_dg_find_peaks_new(const struct gkyl_array_dg_find_peaks_inp *find_peaks_inp,
  const struct gkyl_array *in)
{
  struct gkyl_array_dg_find_peaks *up = gkyl_malloc(sizeof(*up));

  // Copy input parameters.
  up->grid = *find_peaks_inp->grid;
  up->basis = *find_peaks_inp->basis;
  up->range = *find_peaks_inp->range;
  up->range_ext = *find_peaks_inp->range_ext;
  up->search_dir = find_peaks_inp->search_dir;
  up->use_gpu = find_peaks_inp->use_gpu;

  int ndim = find_peaks_inp->grid->ndim;
  int poly_order = find_peaks_inp->basis->poly_order;
  int out_dim = ndim - 1;

  assert(find_peaks_inp->search_dir >= 0 && find_peaks_inp->search_dir < ndim);
  assert(poly_order == 1); // gkyl_array_dg_find_peaks: only p=1 is supported

  // Set up output grid/basis/range.
  if (out_dim == 0) {
    // 1D -> 0D case.
    int cells_1d[1] = { 1 };
    double lower_1d[1] = { 0.0 };
    double upper_1d[1] = { 1.0 };
    gkyl_rect_grid_init(&up->out_grid, 1, lower_1d, upper_1d, cells_1d);
    gkyl_range_init(&up->out_range, 1, (int[]){ 1 }, (int[]){ 1 });
    gkyl_range_init(&up->out_range_ext, 1, (int[]){ 0 }, (int[]){ 2 });
    gkyl_cart_modal_serendip(&up->out_basis, 1, 0);

    int nodes_shape[1] = { 1 };
    gkyl_range_init_from_shape(&up->out_nrange, 1, nodes_shape);
  }
  else if (out_dim == 1) {
    // 2D -> 1D case.
    int preserved_dir = (find_peaks_inp->search_dir == 0) ? 1 : 0;

    int cells_out = find_peaks_inp->grid->cells[preserved_dir];
    double lower_out = find_peaks_inp->grid->lower[preserved_dir];
    double upper_out = find_peaks_inp->grid->upper[preserved_dir];

    gkyl_rect_grid_init(&up->out_grid, 1, &lower_out, &upper_out, &cells_out);

    int lower_idx[1] = { find_peaks_inp->range->lower[preserved_dir] };
    int upper_idx[1] = { find_peaks_inp->range->upper[preserved_dir] };
    gkyl_range_init(&up->out_range, 1, lower_idx, upper_idx);

    int lower_ext_idx[1] = { find_peaks_inp->range_ext->lower[preserved_dir] };
    int upper_ext_idx[1] = { find_peaks_inp->range_ext->upper[preserved_dir] };
    gkyl_range_init(&up->out_range_ext, 1, lower_ext_idx, upper_ext_idx);

    gkyl_cart_modal_serendip(&up->out_basis, 1, poly_order);

    int num_nodes = gkyl_range_shape(&up->out_range, 0) + 1;
    int nodes_shape[1] = {num_nodes};
    gkyl_range_init_from_shape(&up->out_nrange, 1, nodes_shape);
  }
  else {
    assert(false); // dg_find_peaks: only 1D->0D and 2D->1D supported
  }

  // Store node locations for input basis.
  up->nodes = gkyl_array_new(GKYL_DOUBLE, ndim, find_peaks_inp->basis->num_basis);
  find_peaks_inp->basis->node_list(gkyl_array_fetch(up->nodes, 0));

  // Create nodal-to-modal converter.
  up->n2m = gkyl_nodal_ops_new(&up->out_basis, &up->out_grid, false);

  // No device basis on CPU.
  up->out_basis_on_dev = NULL;

  // Compute total_nodes_search for the struct.
  int num_cells_search = find_peaks_inp->range->upper[find_peaks_inp->search_dir]
    - find_peaks_inp->range->lower[find_peaks_inp->search_dir] + 1;
  up->total_nodes_search = num_cells_search + 1;

  // Pre-allocate search-direction working buffers (reused by advance).
  up->search_vals = gkyl_malloc(sizeof(double) * up->total_nodes_search);
  up->search_coords = gkyl_malloc(sizeof(double) * up->total_nodes_search);
  up->search_visited = gkyl_malloc(sizeof(bool) * up->total_nodes_search);

  // Count peaks at middle preserved coordinate.
  int mid_preserved_idx = 0;
  if (out_dim == 1) {
    int preserved_dir = (find_peaks_inp->search_dir == 0) ? 1 : 0;
    mid_preserved_idx = (find_peaks_inp->range->lower[preserved_dir] +
      find_peaks_inp->range->upper[preserved_dir]) / 2;
  }

  // Copy input to host if needed.
  if (up->use_gpu) {
    struct gkyl_array *field_ho = gkyl_array_new(GKYL_DOUBLE, in->ncomp, in->size);
    gkyl_array_copy(field_ho, in);
    count_peaks_along_dir(up, field_ho, mid_preserved_idx, &up->num_peaks, up->peak_types);
    gkyl_array_release(field_ho);
  }
  else {
    count_peaks_along_dir(up, in, mid_preserved_idx, &up->num_peaks, up->peak_types);
  }

  // Allocate output arrays for each peak.
  for (int p = 0; p < up->num_peaks; p++) {
    up->out_vals[p] = gkyl_array_new(GKYL_DOUBLE, up->out_basis.num_basis,
      up->out_range_ext.volume);
    up->out_coords[p] = gkyl_array_new(GKYL_DOUBLE, up->out_basis.num_basis,
      up->out_range_ext.volume);
    up->out_vals_nodal[p] = gkyl_array_new(GKYL_DOUBLE, 1, up->out_nrange.volume);
    up->out_coords_nodal[p] = gkyl_array_new(GKYL_DOUBLE, 1, up->out_nrange.volume);
    up->out_eval_at_peaks_vals_nodal[p] = gkyl_array_new(GKYL_DOUBLE, 1, up->out_nrange.volume);
  }

  // Initialize unused peak arrays to NULL.
  for (int p = up->num_peaks; p < GKYL_DG_FIND_PEAKS_MAX; p++) {
    up->out_vals[p] = NULL;
    up->out_coords[p] = NULL;
    up->out_vals_nodal[p] = NULL;
    up->out_coords_nodal[p] = NULL;
    up->out_eval_at_peaks_vals_nodal[p] = NULL;
  }

  up->flags = 0;
  GKYL_CLEAR_CU_ALLOC(up->flags);
  up->ref_count = gkyl_ref_count_init(gkyl_array_dg_find_peaks_free);
  up->on_dev = up; // CPU object points to itself.

  struct gkyl_array_dg_find_peaks *up_out = up;
#ifdef GKYL_HAVE_CUDA
  if (up->use_gpu) {
    up_out = gkyl_array_dg_find_peaks_new_cu(up);
    gkyl_array_dg_find_peaks_release(up);
  }
#endif

  return up_out;
}

void
gkyl_array_dg_find_peaks_advance(struct gkyl_array_dg_find_peaks *up, const struct gkyl_array *in)
{
#ifdef GKYL_HAVE_CUDA
  if (up->use_gpu) {
    gkyl_array_dg_find_peaks_advance_cu(up, in);
    return;
  }
#endif

  int ndim = up->grid.ndim;
  int out_dim = ndim - 1;

  // Find peaks for each preserved-direction node.
  int num_nodes_out = up->out_nrange.volume;
  for (int pres_node = 0; pres_node < num_nodes_out; pres_node++) {
    find_peaks_for_preserved_node(up, in, pres_node);
  }

  // Transform nodal to modal for each peak.
  if (out_dim == 0) {
    // 1D -> 0D case: modal = nodal (p=0 has no nodal_to_modal function).
    for (int p = 0; p < up->num_peaks; p++) {
      double *val_m = gkyl_array_fetch(up->out_vals[p], 0);
      double *coord_m = gkyl_array_fetch(up->out_coords[p], 0);
      const double *val_n = gkyl_array_cfetch(up->out_vals_nodal[p], 0);
      const double *coord_n = gkyl_array_cfetch(up->out_coords_nodal[p], 0);
      val_m[0] = val_n[0];
      coord_m[0] = coord_n[0];
    }
  }
  else {
    // 2D -> 1D case: use nodal-to-modal transform.
    for (int p = 0; p < up->num_peaks; p++) {
      gkyl_nodal_ops_n2m(up->n2m, &up->out_basis, &up->out_grid,
        &up->out_nrange, &up->out_range, 1, up->out_vals_nodal[p], up->out_vals[p], false);
      gkyl_nodal_ops_n2m(up->n2m, &up->out_basis, &up->out_grid,
        &up->out_nrange, &up->out_range, 1, up->out_coords_nodal[p], up->out_coords[p], false);
    }
  }
}

int
gkyl_array_dg_find_peaks_num_peaks(const struct gkyl_array_dg_find_peaks *up)
{
  return up->num_peaks;
}

enum gkyl_peak_type
gkyl_array_dg_find_peaks_get_type(const struct gkyl_array_dg_find_peaks *up, int peak_idx)
{
  assert(peak_idx >= 0 && peak_idx < up->num_peaks);
  return up->peak_types[peak_idx];
}

const struct gkyl_basis*
gkyl_array_dg_find_peaks_get_basis(const struct gkyl_array_dg_find_peaks *up)
{
  return &up->out_basis;
}

const struct gkyl_rect_grid*
gkyl_array_dg_find_peaks_get_grid(const struct gkyl_array_dg_find_peaks *up)
{
  return &up->out_grid;
}

const struct gkyl_range*
gkyl_array_dg_find_peaks_get_range(const struct gkyl_array_dg_find_peaks *up)
{
  return &up->out_range;
}

const struct gkyl_range*
gkyl_array_dg_find_peaks_get_range_ext(const struct gkyl_array_dg_find_peaks *up)
{
  return &up->out_range_ext;
}

const struct gkyl_range*
gkyl_array_dg_find_peaks_get_nodal_range(const struct gkyl_array_dg_find_peaks *up)
{
  return &up->out_nrange;
}

const struct gkyl_array*
gkyl_array_dg_find_peaks_acquire_vals(const struct gkyl_array_dg_find_peaks *up, int peak_idx)
{
  assert(peak_idx >= 0 && peak_idx < up->num_peaks);
  return gkyl_array_acquire(up->out_vals[peak_idx]);
}

const struct gkyl_array*
gkyl_array_dg_find_peaks_acquire_vals_nodal(const struct gkyl_array_dg_find_peaks *up, int peak_idx)
{
  assert(peak_idx >= 0 && peak_idx < up->num_peaks);
  return gkyl_array_acquire(up->out_vals_nodal[peak_idx]);
}

const struct gkyl_array*
gkyl_array_dg_find_peaks_acquire_coords(const struct gkyl_array_dg_find_peaks *up, int peak_idx)
{
  assert(peak_idx >= 0 && peak_idx < up->num_peaks);
  return gkyl_array_acquire(up->out_coords[peak_idx]);
}

const struct gkyl_array*
gkyl_array_dg_find_peaks_acquire_coords_nodal(const struct gkyl_array_dg_find_peaks *up,
  int peak_idx)
{
  assert(peak_idx >= 0 && peak_idx < up->num_peaks);
  return gkyl_array_acquire(up->out_coords_nodal[peak_idx]);
}

void
gkyl_array_dg_find_peaks_project_on_peaks(struct gkyl_array_dg_find_peaks *up,
  const struct gkyl_array *in_array, struct gkyl_array **out_vals)
{
#ifdef GKYL_HAVE_CUDA
  if (up->use_gpu) {
    gkyl_array_dg_find_peaks_project_on_peaks_cu(up, in_array, out_vals);
    return;
  }
#endif

  int ndim = up->grid.ndim;
  int out_dim = ndim - 1;

  // Evaluate the input array at peak locations for each preserved-direction node.
  int num_nodes_out = up->out_nrange.volume;
  for (int pres_node = 0; pres_node < num_nodes_out; pres_node++) {
    for (int p = 0; p < up->num_peaks; p++) {
      eval_array_at_peaks_for_preserved_node(up, in_array, pres_node,
        up->out_eval_at_peaks_vals_nodal, p);
    }
  }
  // Transform nodal to modal for each peak.
  if (out_dim == 0) {
    // 1D -> 0D case: modal = nodal (p=0 has no nodal_to_modal function).
    for (int p = 0; p < up->num_peaks; p++) {
      double *val_m = gkyl_array_fetch(out_vals[p], 0);
      const double *val_n = gkyl_array_cfetch(up->out_eval_at_peaks_vals_nodal[p], 0);
      val_m[0] = val_n[0];
    }
  }
  else {
    // 2D -> 1D case: use nodal-to-modal transform.
    for (int p = 0; p < up->num_peaks; p++) {
      gkyl_nodal_ops_n2m(up->n2m, &up->out_basis, &up->out_grid,
        &up->out_nrange, &up->out_range, 1, up->out_eval_at_peaks_vals_nodal[p], out_vals[p],
        false);
    }
  }
}

void
gkyl_array_dg_find_peaks_project_on_peak_idx(struct gkyl_array_dg_find_peaks *up,
  const struct gkyl_array *in_array, int peak_idx, struct gkyl_array *out_val)
{
#ifdef GKYL_HAVE_CUDA
  if (up->use_gpu) {
    gkyl_array_dg_find_peaks_project_on_peak_idx_cu(up, in_array, peak_idx, out_val);
    return;
  }
#endif

  int ndim = up->grid.ndim;
  int out_dim = ndim - 1;

  // Evaluate the input array at peak locations for each preserved-direction node.
  int num_nodes_out = up->out_nrange.volume;

  for (int pres_node = 0; pres_node < num_nodes_out; pres_node++) {
    eval_array_at_peaks_for_preserved_node(up, in_array, pres_node,
      up->out_eval_at_peaks_vals_nodal, peak_idx);
  }

  // Transform nodal to modal for each peak.
  if (out_dim == 0) {
    // 1D -> 0D case: modal = nodal (p=0 has no nodal_to_modal function).
    double *val_m = gkyl_array_fetch(out_val, 0);
    const double *val_n = gkyl_array_cfetch(up->out_eval_at_peaks_vals_nodal[peak_idx], 0);
    val_m[0] = val_n[0];
  }
  else {
    // 2D -> 1D case: use nodal-to-modal transform.
    gkyl_nodal_ops_n2m(up->n2m, &up->out_basis, &up->out_grid,
      &up->out_nrange, &up->out_range, 1, up->out_eval_at_peaks_vals_nodal[peak_idx], out_val,
      false);
  }
}

struct gkyl_array_dg_find_peaks*
gkyl_array_dg_find_peaks_acquire(const struct gkyl_array_dg_find_peaks *up)
{
  gkyl_ref_count_inc(&up->ref_count);
  return (struct gkyl_array_dg_find_peaks *)up;
}

void
gkyl_array_dg_find_peaks_free(const struct gkyl_ref_count *ref)
{
  struct gkyl_array_dg_find_peaks *up =
    container_of(ref, struct gkyl_array_dg_find_peaks, ref_count);

  for (int p = 0; p < up->num_peaks; p++) {
    gkyl_array_release(up->out_vals[p]);
    gkyl_array_release(up->out_coords[p]);
    gkyl_array_release(up->out_vals_nodal[p]);
    gkyl_array_release(up->out_coords_nodal[p]);
    gkyl_array_release(up->out_eval_at_peaks_vals_nodal[p]);
  }
  gkyl_array_release(up->nodes);
  gkyl_nodal_ops_release(up->n2m);

  if (GKYL_IS_CU_ALLOC(up->flags)) {
    gkyl_cart_modal_basis_release_cu(up->out_basis_on_dev);
    gkyl_cu_free(up->search_vals);
    gkyl_cu_free(up->search_coords);
    gkyl_cu_free(up->search_visited);
    gkyl_cu_free(up->on_dev);
  }
  else {
    gkyl_free(up->search_vals);
    gkyl_free(up->search_coords);
    gkyl_free(up->search_visited);
  }

  gkyl_free(up);
}

void
gkyl_array_dg_find_peaks_release(struct gkyl_array_dg_find_peaks *up)
{
  gkyl_ref_count_dec(&up->ref_count);
}