ceed-qfunction.c

// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED:  http://github.com/ceed

#include <ceed-impl.h>
#include <ceed.h>
#include <ceed/backend.h>
#include <ceed/jit-tools.h>
#include <limits.h>
#include <stdbool.h>
#include <stdio.h>
#include <string.h>

/// @file
/// Implementation of public CeedQFunction interfaces

/// @cond DOXYGEN_SKIP
static struct CeedQFunction_private ceed_qfunction_none;
/// @endcond

/// @addtogroup CeedQFunctionUser
/// @{

// Indicate that no `CeedQFunction` is provided by the user
const CeedQFunction CEED_QFUNCTION_NONE = &ceed_qfunction_none;

/// @}

/// @cond DOXYGEN_SKIP
static struct {
  char              name[CEED_MAX_RESOURCE_LEN];
  char              source[CEED_MAX_RESOURCE_LEN];
  CeedInt           vec_length;
  CeedQFunctionUser f;
  int (*init)(Ceed ceed, const char *name, CeedQFunction qf);
} gallery_qfunctions[1024];
static size_t num_qfunctions;
/// @endcond

/// ----------------------------------------------------------------------------
/// CeedQFunction Library Internal Functions
/// ----------------------------------------------------------------------------
/// @addtogroup CeedQFunctionDeveloper
/// @{

/**
  @brief Register a gallery `CeedQFunction`

  @param[in] name       Name for this backend to respond to
  @param[in] source     Absolute path to source of `CeedQFunction`, "\path\CEED_DIR\gallery\folder\file.h:function_name"
  @param[in] vec_length Vector length.
                          Caller must ensure that number of quadrature points is a multiple of `vec_length`.
  @param[in] f          Function pointer to evaluate action at quadrature points.
                          See `CeedQFunctionUser`.
  @param[in] init       Initialization function called by @ref CeedQFunctionCreateInteriorByName() when the `CeedQFunction` is selected.

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
int CeedQFunctionRegister(const char *name, const char *source, CeedInt vec_length, CeedQFunctionUser f,
                          int (*init)(Ceed, const char *, CeedQFunction)) {
  const char *relative_file_path;
  int         ierr = 0;

  CeedDebugEnv("Gallery Register: %s", name);
  CeedCall(CeedGetJitRelativePath(source, &relative_file_path));
  CeedPragmaCritical(CeedQFunctionRegister) {
    if (num_qfunctions < sizeof(gallery_qfunctions) / sizeof(gallery_qfunctions[0])) {
      strncpy(gallery_qfunctions[num_qfunctions].name, name, CEED_MAX_RESOURCE_LEN);
      gallery_qfunctions[num_qfunctions].name[CEED_MAX_RESOURCE_LEN - 1] = 0;
      strncpy(gallery_qfunctions[num_qfunctions].source, relative_file_path, CEED_MAX_RESOURCE_LEN);
      gallery_qfunctions[num_qfunctions].source[CEED_MAX_RESOURCE_LEN - 1] = 0;
      gallery_qfunctions[num_qfunctions].vec_length                        = vec_length;
      gallery_qfunctions[num_qfunctions].f                                 = f;
      gallery_qfunctions[num_qfunctions].init                              = init;
      num_qfunctions++;
    } else {
      ierr = 1;
    }
  }
  CeedCheck(ierr == 0, NULL, CEED_ERROR_MAJOR, "Too many gallery CeedQFunctions");
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set a `CeedQFunction` field, used by @ref CeedQFunctionAddInput() and @ref CeedQFunctionAddOutput()

  @param[out] f           `CeedQFunctionField`
  @param[in]  field_name  Name of `CeedQFunction` field
  @param[in]  size        Size of `CeedQFunction` field, (`num_comp * 1`) for @ref CEED_EVAL_NONE and @ref CEED_EVAL_WEIGHT, (`num_comp * 1`) for @ref CEED_EVAL_INTERP for an \f$H^1\f$ space or (`num_comp * dim`) for an \f$H(\mathrm{div})\f$ or \f$H(\mathrm{curl})\f$ space, (`num_comp * dim`) for @ref CEED_EVAL_GRAD, or (num_comp * 1) for @ref CEED_EVAL_DIV, and (`num_comp * curl_dim`) with `curl_dim = 1` if `dim < 3` and `curl_dim = dim` for @ref CEED_EVAL_CURL.
  @param[in]  eval_mode   @ref CEED_EVAL_NONE to use values directly,
                            @ref CEED_EVAL_WEIGHT to use quadrature weights,
                            @ref CEED_EVAL_INTERP to use interpolated values,
                            @ref CEED_EVAL_GRAD to use gradients,
                            @ref CEED_EVAL_DIV to use divergence,
                            @ref CEED_EVAL_CURL to use curl

  @return An error code: 0 - success, otherwise - failure

  @ref Developer
**/
static int CeedQFunctionFieldSet(CeedQFunctionField *f, const char *field_name, CeedInt size, CeedEvalMode eval_mode) {
  CeedCall(CeedCalloc(1, f));
  CeedCall(CeedStringAllocCopy(field_name, (char **)&(*f)->field_name));
  (*f)->size      = size;
  (*f)->eval_mode = eval_mode;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a field of a `CeedQFunction`

  @param[in] field        `CeedQFunction` field to view
  @param[in] field_number Number of field being viewed
  @param[in] in           true for input field, false for output
  @param[in] stream       Stream to view to, e.g., `stdout`

  @return An error code: 0 - success, otherwise - failure

  @ref Utility
**/
static int CeedQFunctionFieldView(CeedQFunctionField field, CeedInt field_number, bool in, FILE *stream) {
  const char  *inout = in ? "Input" : "Output";
  const char  *field_name;
  CeedInt      size;
  CeedEvalMode eval_mode;

  CeedCall(CeedQFunctionFieldGetData(field, &field_name, &size, &eval_mode));
  fprintf(stream,
          "    %s field %" CeedInt_FMT
          ":\n"
          "      Name: \"%s\"\n"
          "      Size: %" CeedInt_FMT
          "\n"
          "      EvalMode: \"%s\"\n",
          inout, field_number, field_name, size, CeedEvalModes[eval_mode]);
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set flag to determine if Fortran interface is used

  @param[in,out] qf     CeedQFunction
  @param[in]     status Boolean value to set as Fortran status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionSetFortranStatus(CeedQFunction qf, bool status) {
  qf->is_fortran = status;
  return CEED_ERROR_SUCCESS;
}

/// @}

/// ----------------------------------------------------------------------------
/// CeedQFunction Backend API
/// ----------------------------------------------------------------------------
/// @addtogroup CeedQFunctionBackend
/// @{

/**
  @brief Get the vector length of a `CeedQFunction`

  @param[in]  qf         `CeedQFunction`
  @param[out] vec_length Variable to store vector length

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetVectorLength(CeedQFunction qf, CeedInt *vec_length) {
  *vec_length = qf->vec_length;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the number of inputs and outputs to a `CeedQFunction`

  @param[in]  qf         `CeedQFunction`
  @param[out] num_input  Variable to store number of input fields
  @param[out] num_output Variable to store number of output fields

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetNumArgs(CeedQFunction qf, CeedInt *num_input, CeedInt *num_output) {
  if (num_input) *num_input = qf->num_input_fields;
  if (num_output) *num_output = qf->num_output_fields;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the name of the user function for a `CeedQFunction`

  @param[in]  qf          `CeedQFunction`
  @param[out] kernel_name Variable to store source path string

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetKernelName(CeedQFunction qf, const char **kernel_name) {
  if (!qf->kernel_name) {
    char *kernel_name_copy;

    if (qf->user_source) {
      const char *kernel_name     = strrchr(qf->user_source, ':') + 1;
      size_t      kernel_name_len = strlen(kernel_name);

      CeedCall(CeedCalloc(kernel_name_len + 1, &kernel_name_copy));
      memcpy(kernel_name_copy, kernel_name, kernel_name_len);
    } else {
      CeedCall(CeedCalloc(1, &kernel_name_copy));
    }
    qf->kernel_name = kernel_name_copy;
  }

  *kernel_name = qf->kernel_name;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the source path string for a `CeedQFunction`

  @param[in]  qf          `CeedQFunction`
  @param[out] source_path Variable to store source path string

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetSourcePath(CeedQFunction qf, const char **source_path) {
  if (!qf->source_path && qf->user_source) {
    Ceed        ceed;
    bool        is_absolute_path;
    char       *source_path_copy;
    const char *absolute_path;
    const char *kernel_name     = strrchr(qf->user_source, ':') + 1;
    size_t      kernel_name_len = strlen(kernel_name);

    CeedCall(CeedQFunctionGetCeed(qf, &ceed));
    CeedCall(CeedCheckFilePath(ceed, qf->user_source, &is_absolute_path));
    if (is_absolute_path) {
      absolute_path = (char *)qf->user_source;
    } else {
      CeedCall(CeedGetJitAbsolutePath(ceed, qf->user_source, &absolute_path));
    }

    size_t source_len = strlen(absolute_path) - kernel_name_len - 1;

    CeedCall(CeedCalloc(source_len + 1, &source_path_copy));
    memcpy(source_path_copy, absolute_path, source_len);
    qf->source_path = source_path_copy;

    if (!is_absolute_path) CeedCall(CeedFree(&absolute_path));
  }

  *source_path = (char *)qf->source_path;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Initialize and load `CeedQFunction` source file into string buffer, including full text of local files in place of `#include "local.h"`.

  The `buffer` is set to `NULL` if there is no `CeedQFunction` source file.

  Note: This function may as well return a mutable buffer, but all current uses
  do not modify it. (This is just a downside of `const` semantics with output
  arguments instead of returns.)

  Note: Caller is responsible for freeing the string buffer with @ref CeedFree().

  @param[in]  qf            `CeedQFunction`
  @param[out] source_buffer String buffer for source file contents

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionLoadSourceToBuffer(CeedQFunction qf, const char **source_buffer) {
  const char *source_path;

  CeedCall(CeedQFunctionGetSourcePath(qf, &source_path));
  *source_buffer = NULL;
  if (source_path) {
    Ceed  ceed;
    char *buffer = NULL;

    CeedCall(CeedQFunctionGetCeed(qf, &ceed));
    CeedCall(CeedLoadSourceToBuffer(ceed, source_path, &buffer));
    *source_buffer = buffer;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the User Function for a `CeedQFunction`

  @param[in]  qf `CeedQFunction`
  @param[out] f  Variable to store user function

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetUserFunction(CeedQFunction qf, CeedQFunctionUser *f) {
  *f = qf->function;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get global context for a `CeedQFunction`.

  Note: For `CeedQFunction` from the Fortran interface, this function will return the Fortran context `CeedQFunctionContext`.

  @param[in]  qf  CeedQFunction
  @param[out] ctx Variable to store CeedQFunctionContext

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetContext(CeedQFunction qf, CeedQFunctionContext *ctx) {
  *ctx = qf->ctx;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get context data of a `CeedQFunction`

  @param[in]  qf       `CeedQFunction`
  @param[in]  mem_type Memory type on which to access the data.
                         If the backend uses a different memory type, this will perform a copy.
  @param[out] data     Data on memory type mem_type

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetContextData(CeedQFunction qf, CeedMemType mem_type, void *data) {
  bool                 is_writable;
  CeedQFunctionContext ctx;

  CeedCall(CeedQFunctionGetContext(qf, &ctx));
  if (ctx) {
    CeedCall(CeedQFunctionIsContextWritable(qf, &is_writable));
    if (is_writable) {
      CeedCall(CeedQFunctionContextGetData(ctx, mem_type, data));
    } else {
      CeedCall(CeedQFunctionContextGetDataRead(ctx, mem_type, data));
    }
  } else {
    *(void **)data = NULL;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Restore context data of a `CeedQFunction`

  @param[in]     qf   `CeedQFunction`
  @param[in,out] data Data to restore

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionRestoreContextData(CeedQFunction qf, void *data) {
  bool                 is_writable;
  CeedQFunctionContext ctx;

  CeedCall(CeedQFunctionGetContext(qf, &ctx));
  if (ctx) {
    CeedCall(CeedQFunctionIsContextWritable(qf, &is_writable));
    if (is_writable) {
      CeedCall(CeedQFunctionContextRestoreData(ctx, data));
    } else {
      CeedCall(CeedQFunctionContextRestoreDataRead(ctx, data));
    }
  }
  *(void **)data = NULL;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get true user context for a `CeedQFunction`

  Note: For all `CeedQFunction` this function will return the user `CeedQFunctionContext` and not interface context `CeedQFunctionContext`, if any such object exists.

  @param[in]  qf  `CeedQFunction`
  @param[out] ctx Variable to store `CeedQFunctionContext`

  @return An error code: 0 - success, otherwise - failure
  @ref Backend
**/
int CeedQFunctionGetInnerContext(CeedQFunction qf, CeedQFunctionContext *ctx) {
  CeedQFunctionContext qf_ctx;

  CeedCall(CeedQFunctionGetContext(qf, &qf_ctx));
  if (qf->is_fortran) {
    CeedFortranContext fortran_ctx = NULL;

    CeedCall(CeedQFunctionContextGetData(qf_ctx, CEED_MEM_HOST, &fortran_ctx));
    *ctx = fortran_ctx->inner_ctx;
    CeedCall(CeedQFunctionContextRestoreData(qf_ctx, &fortran_ctx));
  } else {
    *ctx = qf_ctx;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get inner context data of a `CeedQFunction`

  @param[in]  qf       `CeedQFunction`
  @param[in]  mem_type Memory type on which to access the data.
                         If the backend uses a different memory type, this will perform a copy.
  @param[out] data     Data on memory type mem_type

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetInnerContextData(CeedQFunction qf, CeedMemType mem_type, void *data) {
  bool                 is_writable;
  CeedQFunctionContext ctx;

  CeedCall(CeedQFunctionGetInnerContext(qf, &ctx));
  if (ctx) {
    CeedCall(CeedQFunctionIsContextWritable(qf, &is_writable));
    if (is_writable) {
      CeedCall(CeedQFunctionContextGetData(ctx, mem_type, data));
    } else {
      CeedCall(CeedQFunctionContextGetDataRead(ctx, mem_type, data));
    }
  } else {
    *(void **)data = NULL;
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Restore inner context data of a `CeedQFunction`

  @param[in]     qf   `CeedQFunction`
  @param[in,out] data Data to restore

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionRestoreInnerContextData(CeedQFunction qf, void *data) {
  bool                 is_writable;
  CeedQFunctionContext ctx;

  CeedCall(CeedQFunctionGetInnerContext(qf, &ctx));
  if (ctx) {
    CeedCall(CeedQFunctionIsContextWritable(qf, &is_writable));
    if (is_writable) {
      CeedCall(CeedQFunctionContextRestoreData(ctx, data));
    } else {
      CeedCall(CeedQFunctionContextRestoreDataRead(ctx, data));
    }
  }
  *(void **)data = NULL;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Determine if `CeedQFunction` is identity

  @param[in]  qf          `CeedQFunction`
  @param[out] is_identity Variable to store identity status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionIsIdentity(CeedQFunction qf, bool *is_identity) {
  *is_identity = qf->is_identity;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Determine if `CeedQFunctionContext` is writable

  @param[in]  qf          `CeedQFunction`
  @param[out] is_writable Variable to store context writeable status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionIsContextWritable(CeedQFunction qf, bool *is_writable) {
  *is_writable = qf->is_context_writable;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get backend data of a `CeedQFunction`

  @param[in]  qf   `CeedQFunction`
  @param[out] data Variable to store data

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionGetData(CeedQFunction qf, void *data) {
  *(void **)data = qf->data;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set backend data of a `CeedQFunction`

  @param[in,out] qf   `CeedQFunction`
  @param[in]     data Data to set

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionSetData(CeedQFunction qf, void *data) {
  qf->data = data;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get a boolean value indicating if the `CeedQFunction` is immutable

  @param[in]  qf           `CeedOperator`
  @param[out] is_immutable Variable to store immutability status

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionIsImmutable(CeedQFunction qf, bool *is_immutable) {
  *is_immutable = qf->is_immutable;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set the immutable flag of a `CeedQFunction` to `true`

  @param[in,out] qf `CeedQFunction`

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionSetImmutable(CeedQFunction qf) {
  qf->is_immutable = true;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Increment the reference counter for a `CeedQFunction`

  @param[in,out] qf `CeedQFunction` to increment the reference counter

  @return An error code: 0 - success, otherwise - failure

  @ref Backend
**/
int CeedQFunctionReference(CeedQFunction qf) {
  qf->ref_count++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Estimate number of FLOPs per quadrature required to apply `CeedQFunction`

  @param[in]  qf    `CeedQFunction` to estimate FLOPs for
  @param[out] flops Address of variable to hold FLOPs estimate

  @ref Backend
**/
int CeedQFunctionGetFlopsEstimate(CeedQFunction qf, CeedSize *flops) {
  *flops = qf->user_flop_estimate;
  return CEED_ERROR_SUCCESS;
}

/// @}

/// ----------------------------------------------------------------------------
/// CeedQFunction Public API
/// ----------------------------------------------------------------------------
/// @addtogroup CeedQFunctionUser
/// @{

/**
  @brief Create a `CeedQFunction` for evaluating interior (volumetric) terms

  @param[in]  ceed       `Ceed` object used to create the `CeedQFunction`
  @param[in]  vec_length Vector length.
                           Caller must ensure that number of quadrature points is a multiple of `vec_length`.
  @param[in]  f          Function pointer to evaluate action at quadrature points.
                           See `CeedQFunctionUser`.
  @param[in]  source     Absolute path to source of `CeedQFunctionUser`, "\abs_path\file.h:function_name".
                           The entire source file must only contain constructs supported by all targeted backends (i.e. CUDA for `/gpu/cuda`, OpenCL/SYCL for `/gpu/sycl`, etc.).
                           The entire contents of this file and all locally included files are used during JiT compilation for GPU backends.
                           All source files must be at the provided filepath at runtime for JiT to function.
  @param[out] qf         Address of the variable where the newly created `CeedQFunction` will be stored

  @return An error code: 0 - success, otherwise - failure

  See \ref CeedQFunctionUser for details on the call-back function `f` arguments.

  @ref User
**/
int CeedQFunctionCreateInterior(Ceed ceed, CeedInt vec_length, CeedQFunctionUser f, const char *source, CeedQFunction *qf) {
  char *user_source_copy;

  if (!ceed->QFunctionCreate) {
    Ceed delegate;

    CeedCall(CeedGetObjectDelegate(ceed, &delegate, "QFunction"));
    CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedQFunctionCreateInterior");
    CeedCall(CeedQFunctionCreateInterior(delegate, vec_length, f, source, qf));
    return CEED_ERROR_SUCCESS;
  }

  CeedCheck(!strlen(source) || strrchr(source, ':'), ceed, CEED_ERROR_INCOMPLETE,
            "Provided path to source does not include function name. Provided: \"%s\"\nRequired: \"\\abs_path\\file.h:function_name\"", source);

  CeedCall(CeedCalloc(1, qf));
  CeedCall(CeedReferenceCopy(ceed, &(*qf)->ceed));
  (*qf)->ref_count           = 1;
  (*qf)->vec_length          = vec_length;
  (*qf)->is_identity         = false;
  (*qf)->is_context_writable = true;
  (*qf)->function            = f;
  (*qf)->user_flop_estimate  = -1;
  if (strlen(source)) {
    size_t user_source_len = strlen(source);

    CeedCall(CeedCalloc(user_source_len + 1, &user_source_copy));
    memcpy(user_source_copy, source, user_source_len);
    (*qf)->user_source = user_source_copy;
  }
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*qf)->input_fields));
  CeedCall(CeedCalloc(CEED_FIELD_MAX, &(*qf)->output_fields));
  CeedCall(ceed->QFunctionCreate(*qf));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create a `CeedQFunction` for evaluating interior (volumetric) terms by name

  @param[in]  ceed `Ceed` object used to create the `CeedQFunction`
  @param[in]  name Name of `CeedQFunction` to use from gallery
  @param[out] qf   Address of the variable where the newly created `CeedQFunction` will be stored

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionCreateInteriorByName(Ceed ceed, const char *name, CeedQFunction *qf) {
  size_t match_len = 0, match_index = UINT_MAX;

  CeedCall(CeedQFunctionRegisterAll());
  // Find matching backend
  CeedCheck(name, ceed, CEED_ERROR_INCOMPLETE, "No CeedQFunction name provided");
  for (size_t i = 0; i < num_qfunctions; i++) {
    size_t      n;
    const char *curr_name = gallery_qfunctions[i].name;
    for (n = 0; curr_name[n] && curr_name[n] == name[n]; n++) {
    }
    if (n > match_len) {
      match_len   = n;
      match_index = i;
    }
  }
  CeedCheck(match_len > 0, ceed, CEED_ERROR_UNSUPPORTED, "No suitable gallery CeedQFunction");

  // Create QFunction
  CeedCall(CeedQFunctionCreateInterior(ceed, gallery_qfunctions[match_index].vec_length, gallery_qfunctions[match_index].f,
                                       gallery_qfunctions[match_index].source, qf));

  // QFunction specific setup
  CeedCall(gallery_qfunctions[match_index].init(ceed, name, *qf));

  // Copy name
  CeedCall(CeedStringAllocCopy(name, (char **)&(*qf)->gallery_name));
  (*qf)->is_gallery = true;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Create an identity `CeedQFunction`.

  Inputs are written into outputs in the order given.
  This is useful for `CeedOperator that can be represented with only the action of a `CeedElemRestriction` and `CeedBasis`, such as restriction and prolongation operators for p-multigrid.
  Backends may optimize `CeedOperator` with this `CeedQFunction` to avoid the copy of input data to output fields by using the same memory location for both.

  @param[in]  ceed     `Ceed` object used to create the `CeedQFunction`
  @param[in]  size     Size of the `CeedQFunction` fields
  @param[in]  in_mode  @ref CeedEvalMode for input to `CeedQFunction`
  @param[in]  out_mode @ref CeedEvalMode for output to `CeedQFunction`
  @param[out] qf       Address of the variable where the newly created `CeedQFunction` will be stored

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionCreateIdentity(Ceed ceed, CeedInt size, CeedEvalMode in_mode, CeedEvalMode out_mode, CeedQFunction *qf) {
  CeedQFunctionContext  ctx;
  CeedContextFieldLabel size_label;

  CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Identity", qf));
  CeedCall(CeedQFunctionAddInput(*qf, "input", size, in_mode));
  CeedCall(CeedQFunctionAddOutput(*qf, "output", size, out_mode));

  (*qf)->is_identity = true;

  CeedCall(CeedQFunctionGetContext(*qf, &ctx));
  CeedCall(CeedQFunctionContextGetFieldLabel(ctx, "size", &size_label));
  CeedCall(CeedQFunctionContextSetInt32(ctx, size_label, &size));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Copy the pointer to a `CeedQFunction`.

  Both pointers should be destroyed with @ref CeedQFunctionDestroy().

  Note: If the value of `*qf_copy` passed to this function is non-NULL, then it is assumed that `*qf_copy` is a pointer to a `CeedQFunction`.
        This `CeedQFunction` will be destroyed if `*qf_copy` is the only reference to this `CeedQFunction`.

  @param[in]  qf      `CeedQFunction` to copy reference to
  @param[out] qf_copy Variable to store copied reference

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionReferenceCopy(CeedQFunction qf, CeedQFunction *qf_copy) {
  CeedCall(CeedQFunctionReference(qf));
  CeedCall(CeedQFunctionDestroy(qf_copy));
  *qf_copy = qf;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Add a `CeedQFunction` input

  @param[in,out] qf         `CeedQFunction`
  @param[in]     field_name Name of `CeedQFunction` field
  @param[in]     size       Size of `CeedQFunction` field, (`num_comp * 1`) for @ref CEED_EVAL_NONE, (`num_comp * 1`) for @ref CEED_EVAL_INTERP for an \f$H^1\f$ space or (`num_comp * dim`) for an \f$H(\mathrm{div})\f$ or \f$H(\mathrm{curl})\f$ space, (`num_comp * dim`) for @ref CEED_EVAL_GRAD, or (`num_comp * 1`) for @ref CEED_EVAL_DIV, and (`num_comp * curl_dim`) with `curl_dim = 1` if `dim < 3` otherwise `curl_dim = dim` for @ref CEED_EVAL_CURL.
  @param[in]     eval_mode  @ref CEED_EVAL_NONE to use values directly,
                              @ref CEED_EVAL_INTERP to use interpolated values,
                              @ref CEED_EVAL_GRAD to use gradients,
                              @ref CEED_EVAL_DIV to use divergence,
                              @ref CEED_EVAL_CURL to use curl

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionAddInput(CeedQFunction qf, const char *field_name, CeedInt size, CeedEvalMode eval_mode) {
  bool is_immutable;
  Ceed ceed;

  CeedCall(CeedQFunctionGetCeed(qf, &ceed));
  CeedCall(CeedQFunctionIsImmutable(qf, &is_immutable));
  CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "QFunction cannot be changed after set as immutable");
  CeedCheck(eval_mode != CEED_EVAL_WEIGHT || size == 1, ceed, CEED_ERROR_DIMENSION, "CEED_EVAL_WEIGHT should have size 1");
  for (CeedInt i = 0; i < qf->num_input_fields; i++) {
    CeedCheck(strcmp(field_name, qf->input_fields[i]->field_name), ceed, CEED_ERROR_MINOR, "CeedQFunction field names must be unique");
  }
  for (CeedInt i = 0; i < qf->num_output_fields; i++) {
    CeedCheck(strcmp(field_name, qf->output_fields[i]->field_name), ceed, CEED_ERROR_MINOR, "CeedQFunction field names must be unique");
  }
  CeedCall(CeedQFunctionFieldSet(&qf->input_fields[qf->num_input_fields], field_name, size, eval_mode));
  qf->num_input_fields++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Add a `CeedQFunction` output

  @param[in,out] qf         `CeedQFunction`
  @param[in]     field_name Name of `CeedQFunction` field
  @param[in]     size       Size of `CeedQFunction` field, (`num_comp * 1`) for @ref CEED_EVAL_NONE, (`num_comp * 1`) for @ref CEED_EVAL_INTERP for an \f$H^1\f$ space or (`num_comp * dim`) for an \f$H(\mathrm{div})\f$ or \f$H(\mathrm{curl})\f$ space, (`num_comp * dim`) for @ref CEED_EVAL_GRAD, or (`num_comp * 1`) for @ref CEED_EVAL_DIV, and (`num_comp * curl_dim`) with `curl_dim = 1` if `dim < 3` else dim for @ref CEED_EVAL_CURL.
  @param[in]     eval_mode  @ref CEED_EVAL_NONE to use values directly,
                              @ref CEED_EVAL_INTERP to use interpolated values,
                              @ref CEED_EVAL_GRAD to use gradients,
                              @ref CEED_EVAL_DIV to use divergence,
                              @ref CEED_EVAL_CURL to use curl.

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionAddOutput(CeedQFunction qf, const char *field_name, CeedInt size, CeedEvalMode eval_mode) {
  bool is_immutable;
  Ceed ceed;

  CeedCall(CeedQFunctionGetCeed(qf, &ceed));
  CeedCall(CeedQFunctionIsImmutable(qf, &is_immutable));
  CeedCheck(!is_immutable, ceed, CEED_ERROR_MAJOR, "CeedQFunction cannot be changed after set as immutable");
  CeedCheck(eval_mode != CEED_EVAL_WEIGHT, ceed, CEED_ERROR_DIMENSION, "Cannot create CeedQFunction output with CEED_EVAL_WEIGHT");
  for (CeedInt i = 0; i < qf->num_input_fields; i++) {
    CeedCheck(strcmp(field_name, qf->input_fields[i]->field_name), ceed, CEED_ERROR_MINOR, "CeedQFunction field names must be unique");
  }
  for (CeedInt i = 0; i < qf->num_output_fields; i++) {
    CeedCheck(strcmp(field_name, qf->output_fields[i]->field_name), ceed, CEED_ERROR_MINOR, "CeedQFunction field names must be unique");
  }
  CeedCall(CeedQFunctionFieldSet(&qf->output_fields[qf->num_output_fields], field_name, size, eval_mode));
  qf->num_output_fields++;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `CeedQFunctionField` of a `CeedQFunction`

  Note: Calling this function asserts that setup is complete and sets the `CeedQFunction` as immutable.

  @param[in]  qf                `CeedQFunction`
  @param[out] num_input_fields  Variable to store number of input fields
  @param[out] input_fields      Variable to store input fields
  @param[out] num_output_fields Variable to store number of output fields
  @param[out] output_fields     Variable to store output fields

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionGetFields(CeedQFunction qf, CeedInt *num_input_fields, CeedQFunctionField **input_fields, CeedInt *num_output_fields,
                           CeedQFunctionField **output_fields) {
  CeedCall(CeedQFunctionSetImmutable(qf));
  if (num_input_fields) *num_input_fields = qf->num_input_fields;
  if (input_fields) *input_fields = qf->input_fields;
  if (num_output_fields) *num_output_fields = qf->num_output_fields;
  if (output_fields) *output_fields = qf->output_fields;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the name of a `CeedQFunctionField`

  @param[in]  qf_field   `CeedQFunctionField`
  @param[out] field_name Variable to store the field name

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionFieldGetName(CeedQFunctionField qf_field, const char **field_name) {
  *field_name = qf_field->field_name;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the number of components of a `CeedQFunctionField`

  @param[in]  qf_field `CeedQFunctionField`
  @param[out] size     Variable to store the size of the field

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionFieldGetSize(CeedQFunctionField qf_field, CeedInt *size) {
  *size = qf_field->size;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the @ref CeedEvalMode of a `CeedQFunctionField`

  @param[in]  qf_field  `CeedQFunctionField`
  @param[out] eval_mode Variable to store the field evaluation mode

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionFieldGetEvalMode(CeedQFunctionField qf_field, CeedEvalMode *eval_mode) {
  *eval_mode = qf_field->eval_mode;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the data of a `CeedQFunctionField`.

  Any arguments set as `NULL` are ignored.

  @param[in]  qf_field   `CeedQFunctionField`
  @param[out] field_name Variable to store the field name
  @param[out] size       Variable to store the size of the field
  @param[out] eval_mode  Variable to store the field evaluation mode

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionFieldGetData(CeedQFunctionField qf_field, const char **field_name, CeedInt *size, CeedEvalMode *eval_mode) {
  if (field_name) CeedCall(CeedQFunctionFieldGetName(qf_field, field_name));
  if (size) CeedCall(CeedQFunctionFieldGetSize(qf_field, size));
  if (eval_mode) CeedCall(CeedQFunctionFieldGetEvalMode(qf_field, eval_mode));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set global context for a `CeedQFunction`

  @param[in,out] qf  `CeedQFunction`
  @param[in]     ctx Context data to set

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionSetContext(CeedQFunction qf, CeedQFunctionContext ctx) {
  CeedCall(CeedQFunctionContextDestroy(&qf->ctx));
  qf->ctx = ctx;
  if (ctx) CeedCall(CeedQFunctionContextReference(ctx));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set writability of `CeedQFunctionContext` when calling the `CeedQFunctionUser`.

  The default value is `is_writable == true`.

  Setting `is_writable == true` indicates the `CeedQFunctionUser` writes into the `CeedQFunctionContext` and requires memory synchronization after calling @ref CeedQFunctionApply().

  Setting 'is_writable == false' asserts that `CeedQFunctionUser` does not modify the `CeedQFunctionContext`.
  Violating this assertion may lead to inconsistent data.

  Setting `is_writable == false` may offer a performance improvement on GPU backends.

  @param[in,out] qf          `CeedQFunction`
  @param[in]     is_writable Boolean flag for writability status

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionSetContextWritable(CeedQFunction qf, bool is_writable) {
  qf->is_context_writable = is_writable;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Set estimated number of FLOPs per quadrature required to apply `CeedQFunction`

  @param[in]  qf    `CeedQFunction` to estimate FLOPs for
  @param[out] flops FLOPs per quadrature point estimate

  @ref Backend
**/
int CeedQFunctionSetUserFlopsEstimate(CeedQFunction qf, CeedSize flops) {
  CeedCheck(flops >= 0, CeedQFunctionReturnCeed(qf), CEED_ERROR_INCOMPATIBLE, "Must set non-negative FLOPs estimate");
  qf->user_flop_estimate = flops;
  return CEED_ERROR_SUCCESS;
}

/**
  @brief View a `CeedQFunction`

  @param[in] qf     `CeedQFunction` to view
  @param[in] stream Stream to write; typically `stdout` or a file

  @return Error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionView(CeedQFunction qf, FILE *stream) {
  const char *kernel_name;

  CeedCall(CeedQFunctionGetKernelName(qf, &kernel_name));
  fprintf(stream, "%sCeedQFunction - %s\n", qf->is_gallery ? "Gallery " : "User ", qf->is_gallery ? qf->gallery_name : kernel_name);

  fprintf(stream, "  %" CeedInt_FMT " input field%s:\n", qf->num_input_fields, qf->num_input_fields > 1 ? "s" : "");
  for (CeedInt i = 0; i < qf->num_input_fields; i++) {
    CeedCall(CeedQFunctionFieldView(qf->input_fields[i], i, 1, stream));
  }

  fprintf(stream, "  %" CeedInt_FMT " output field%s:\n", qf->num_output_fields, qf->num_output_fields > 1 ? "s" : "");
  for (CeedInt i = 0; i < qf->num_output_fields; i++) {
    CeedCall(CeedQFunctionFieldView(qf->output_fields[i], i, 0, stream));
  }
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Get the `Ceed` associated with a `CeedQFunction`

  @param[in]  qf   `CeedQFunction`
  @param[out] ceed Variable to store`Ceed`

  @return An error code: 0 - success, otherwise - failure

  @ref Advanced
**/
int CeedQFunctionGetCeed(CeedQFunction qf, Ceed *ceed) {
  *ceed = CeedQFunctionReturnCeed(qf);
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Return the `Ceed` associated with a `CeedQFunction`

  @param[in]  qf   `CeedQFunction`

  @return `Ceed` associated with the `qf`

  @ref Advanced
**/
Ceed CeedQFunctionReturnCeed(CeedQFunction qf) { return qf->ceed; }

/**
  @brief Apply the action of a `CeedQFunction`

  Note: Calling this function asserts that setup is complete and sets the `CeedQFunction` as immutable.

  @param[in]  qf `CeedQFunction`
  @param[in]  Q  Number of quadrature points
  @param[in]  u  Array of input `CeedVector`
  @param[out] v  Array of output `CeedVector`

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionApply(CeedQFunction qf, CeedInt Q, CeedVector *u, CeedVector *v) {
  CeedInt vec_length;
  Ceed    ceed;

  CeedCall(CeedQFunctionGetCeed(qf, &ceed));
  CeedCheck(qf->Apply, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedQFunctionApply");
  CeedCall(CeedQFunctionGetVectorLength(qf, &vec_length));
  CeedCheck(Q % vec_length == 0, ceed, CEED_ERROR_DIMENSION, "Number of quadrature points %" CeedInt_FMT " must be a multiple of %" CeedInt_FMT, Q,
            qf->vec_length);
  CeedCall(CeedQFunctionSetImmutable(qf));
  CeedCall(qf->Apply(qf, Q, u, v));
  return CEED_ERROR_SUCCESS;
}

/**
  @brief Destroy a `CeedQFunction`

  @param[in,out] qf `CeedQFunction` to destroy

  @return An error code: 0 - success, otherwise - failure

  @ref User
**/
int CeedQFunctionDestroy(CeedQFunction *qf) {
  if (!*qf || --(*qf)->ref_count > 0) {
    *qf = NULL;
    return CEED_ERROR_SUCCESS;
  }
  // Backend destroy
  if ((*qf)->Destroy) {
    CeedCall((*qf)->Destroy(*qf));
  }
  // Free fields
  for (CeedInt i = 0; i < (*qf)->num_input_fields; i++) {
    CeedCall(CeedFree(&(*(*qf)->input_fields[i]).field_name));
    CeedCall(CeedFree(&(*qf)->input_fields[i]));
  }
  for (CeedInt i = 0; i < (*qf)->num_output_fields; i++) {
    CeedCall(CeedFree(&(*(*qf)->output_fields[i]).field_name));
    CeedCall(CeedFree(&(*qf)->output_fields[i]));
  }
  CeedCall(CeedFree(&(*qf)->input_fields));
  CeedCall(CeedFree(&(*qf)->output_fields));

  // User context data object
  CeedCall(CeedQFunctionContextDestroy(&(*qf)->ctx));

  CeedCall(CeedFree(&(*qf)->user_source));
  CeedCall(CeedFree(&(*qf)->source_path));
  CeedCall(CeedFree(&(*qf)->gallery_name));
  CeedCall(CeedFree(&(*qf)->kernel_name));
  CeedCall(CeedDestroy(&(*qf)->ceed));
  CeedCall(CeedFree(qf));
  return CEED_ERROR_SUCCESS;
}

/// @}