v0.12.0/UserQFunctions.html

<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>Defining User Q-Functions · LibCEED.jl Docs</title><link href="https://fonts.googleapis.com/css?family=Lato|Roboto+Mono" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.0/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.11.1/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="assets/documenter.js"></script><script src="siteinfo.js"></script><script src="../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-dark.css" data-theme-name="documenter-dark"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><div class="docs-package-name"><span class="docs-autofit">LibCEED.jl Docs</span></div><form class="docs-search" action="search.html"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="index.html">Home</a></li><li><span class="tocitem">Ceed Objects</span><ul><li><a class="tocitem" href="Ceed.html">Ceed</a></li><li><a class="tocitem" href="CeedVector.html">CeedVector</a></li><li><a class="tocitem" href="ElemRestriction.html">ElemRestriction</a></li><li><a class="tocitem" href="Basis.html">Basis</a></li><li><a class="tocitem" href="QFunction.html">QFunction</a></li><li><a class="tocitem" href="Operator.html">Operator</a></li></ul></li><li><span class="tocitem">Utilities</span><ul><li><a class="tocitem" href="Misc.html">Linear Algebra</a></li><li><a class="tocitem" href="Globals.html">Constants and Enumerations</a></li><li><a class="tocitem" href="Quadrature.html">Quadrature</a></li></ul></li><li><a class="tocitem" href="LibCEED.html">Library configuration</a></li><li><a class="tocitem" href="C.html">Low-level C interface</a></li><li class="is-active"><a class="tocitem" href="UserQFunctions.html">Defining User Q-Functions</a><ul class="internal"><li><a class="tocitem" href="#Apply-mass-Q-function-in-C"><span>Apply mass Q-function in C</span></a></li><li><a class="tocitem" href="#Apply-mass-Q-function-in-Julia"><span>Apply mass Q-function in Julia</span></a></li><li><a class="tocitem" href="#applydiff"><span>Apply diffusion Q-function in Julia</span></a></li><li><a class="tocitem" href="#GPU-Kernels"><span>GPU Kernels</span></a></li></ul></li><li><a class="tocitem" href="Examples.html">Examples</a></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href="UserQFunctions.html">Defining User Q-Functions</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href="UserQFunctions.html">Defining User Q-Functions</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://github.com/CEED/libCEED/blob/master/julia/LibCEED.jl/docs/src/UserQFunctions.md" title="Edit on GitHub"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitHub</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Defining-User-Q-Functions"><a class="docs-heading-anchor" href="#Defining-User-Q-Functions">Defining User Q-Functions</a><a id="Defining-User-Q-Functions-1"></a><a class="docs-heading-anchor-permalink" href="#Defining-User-Q-Functions" title="Permalink"></a></h1><p>An important feature of LibCEED.jl is the ability to define <a href="https://libceed.org/en/latest/libCEEDapi/#gallery-of-qfunctions">user Q-functions</a> natively in Julia. These user Q-functions work with both the CPU and CUDA backends.</p><p>User Q-functions describe the action of the <span>$D$</span> operator at quadrature points (see <a href="https://libceed.org/en/latest/libCEEDapi/#theoretical-framework">libCEED&#39;s theoretical framework</a>). Since the Q-functions are invoked at every quadrature point, efficiency is very important.</p><h2 id="Apply-mass-Q-function-in-C"><a class="docs-heading-anchor" href="#Apply-mass-Q-function-in-C">Apply mass Q-function in C</a><a id="Apply-mass-Q-function-in-C-1"></a><a class="docs-heading-anchor-permalink" href="#Apply-mass-Q-function-in-C" title="Permalink"></a></h2><p>Before describing how to define user Q-functions in Julia, we will briefly given an example of a user Q-function defined in C. This is the &quot;apply mass&quot; Q-function from <code>ex1-volume.c</code>, which computes the action of the mass operator. The mass operator on each element can be written as <span>$B^\intercal D B$</span>, where <span>$B$</span> is the basis operator, and <span>$D$</span> represents multiplication by quadrature weights and geometric factors (i.e. the determinant of the mesh transformation Jacobian at each qudarture point). It is the action of <span>$D$</span> that the Q-function must implement. The C source of the Q-function is:</p><pre><code class="language-c">/// libCEED Q-function for applying a mass operator
CEED_QFUNCTION(f_apply_mass)(void *ctx, const CeedInt Q,
                             const CeedScalar *const *in,
                             CeedScalar *const *out) {
  const CeedScalar *u = in[0], *qdata = in[1];
  CeedScalar *v = out[0];
  // Quadrature Point Loop
  CeedPragmaSIMD
  for (CeedInt i=0; i&lt;Q; i++) {
    v[i] = qdata[i] * u[i];
  } // End of Quadrature Point Loop
  return 0;
}</code></pre><p>From this example, we see that a user Q-function is a C callback that takes a &quot;data context&quot; pointer, a number of quadrature points, and two arrays of arrays, one for inputs, and one for outputs.</p><p>In this example, the first input array is <code>u</code>, which is the value of the trial function evaluated at each quadrature point. The second input array is <code>qdata</code>, which contains the precomputed geometric factors. There is only one output array, <code>v</code>, which will store the pointwise product of <code>u</code> and <code>data</code>. Given the definition of this Q-function, the <code>CeedQFunction</code> object is created by</p><pre><code class="language-c">CeedQFunctionCreateInterior(ceed, 1, f_apply_mass, f_apply_mass_loc, &amp;apply_qfunc);
CeedQFunctionAddInput(apply_qfunc, &quot;u&quot;, 1, CEED_EVAL_INTERP);
CeedQFunctionAddInput(apply_qfunc, &quot;qdata&quot;, 1, CEED_EVAL_NONE);
CeedQFunctionAddOutput(apply_qfunc, &quot;v&quot;, 1, CEED_EVAL_INTERP);</code></pre><p>When adding the inputs and outputs, <code>CEED_EVAL_INTERP</code> indicates that the <span>$B$</span> basis operator should be used to interpolate the trial and test functions from nodal points to quadrature points, and <code>CEED_EVAL_NONE</code> indicates that the <code>qdata</code> is already precomputed at quadrature points, and no interpolation is requried.</p><h2 id="Apply-mass-Q-function-in-Julia"><a class="docs-heading-anchor" href="#Apply-mass-Q-function-in-Julia">Apply mass Q-function in Julia</a><a id="Apply-mass-Q-function-in-Julia-1"></a><a class="docs-heading-anchor-permalink" href="#Apply-mass-Q-function-in-Julia" title="Permalink"></a></h2><p>We now replicate this Q-function in Julia. The main way of defining user Q-functions in Julia is using the <a href="QFunction.html#LibCEED.@interior_qf"><code>@interior_qf</code></a> macro. The above C code (both the definition of the Q-function, its creation, and adding the inputs and outputs) is analogous to the following Julia code:</p><pre><code class="language-julia">@interior_qf apply_qfunc = (
    ceed, Q,
    (u, :in, EVAL_INTERP, Q), (qdata, :in, EVAL_NONE, Q),
    (v, :out, EVAL_INTERP, Q),
    @inbounds @simd for i=1:Q
        v[i] = qdata[i]*u[i]
    end
)</code></pre><p>This creates a <a href="QFunction.html#LibCEED.QFunction"><code>QFunction</code></a> object named <code>apply_qfunc</code>. The Q-function is defined by the tuple on the right-hand side. <code>ceed</code> is the name of the <a href="Ceed.html#LibCEED.Ceed"><code>Ceed</code></a> object where the Q-function will be created, and the second argument, <code>Q</code>, is the name of that variable that will contain the number of quadrature points. The next three arguments are specifications of the input and output fields:</p><pre><code class="language-julia">    (u, :in, EVAL_INTERP, Q),
    (qdata, :in, EVAL_NONE, Q),
    (v, :out, EVAL_INTERP, Q),</code></pre><p>Each input or output field specification is a tuple, where the first entry is the name of the array, and the second entry is either <code>:in</code> or <code>:out</code>, according to whether the array is an input or output array. The third entry is the <a href="Globals.html#LibCEED.EvalMode"><code>EvalMode</code></a> of the field. The remaining entries are the dimensions of the array. The first dimension is always equal to the number of quadrature points. In this case, all the arrays are simply vectors whose size is equal to the number of quadrature points, but in more sophisticated examples (e.g. the <a href="UserQFunctions.html#applydiff">apply diffusion Q-function</a>) these arrays could consists of vectors or matrices at each quadrature point. After providing all of the array specifications, the body of the Q-function is provided.</p><h2 id="applydiff"><a class="docs-heading-anchor" href="#applydiff">Apply diffusion Q-function in Julia</a><a id="applydiff-1"></a><a class="docs-heading-anchor-permalink" href="#applydiff" title="Permalink"></a></h2><p>For a more sophisticated example of a Q-function, we consider the &quot;apply diffusion&quot; Q-function, used in <code>ex2-surface</code>. This Q-function computes the action of the diffusion operator. When written in the form <span>$B^\intercal D B$</span>, in this case <span>$B$</span> represents the basis gradient matrix, and <span>$D$</span> represents multiplication by <span>$w \det(J) J^{-\intercal} J^{-1}$</span>, where <span>$J$</span> is the mesh transformation Jacobian, and <span>$w$</span> is the quadrature weight.</p><p>This Q-function is implemented in Julia as follows:</p><pre><code class="language-julia">@interior_qf apply_qfunc = (
    ceed, Q, dim=dim,
    (du, :in, EVAL_GRAD, Q, dim),
    (qdata, :in, EVAL_NONE, Q, dim*(dim+1)÷2),
    (dv, :out, EVAL_GRAD, Q, dim),
    @inbounds @simd for i=1:Q
        dXdxdXdxT = getvoigt(@view(qdata[i,:]), CeedDim(dim))
        dui = SVector{dim}(@view(du[i,:]))
        dv[i,:] .= dXdxdXdxT*dui
    end
)</code></pre><p>In contrast to the previous example, before the field specifications, this Q-function includes a <em>constant definition</em> <code>dim=dim</code>. The <a href="QFunction.html#LibCEED.@interior_qf"><code>@interior_qf</code></a> macro allows for any number of constant definitions, which make the specified values available within the body of the Q-function as compile-time constants.</p><p>In this example, <code>dim</code> is either 1, 2, or 3 according to the spatial dimension of the problem. When the user Q-function is defined, LibCEED.jl will JIT compile the body of the Q-function and make it available to libCEED as a C callback. In the body of this Q-function, <code>dim</code> will be available, and its value will be a compile-time constant, allowing for (static) dispatch based on the value of <code>dim</code>, and eliminating branching.</p><p>Note that <code>dim</code> is also available for use in the field specifications. In this example, the field specifications are slightly more involved that in the previous example. The arrays are given by</p><pre><code class="language-julia">    (du, :in, EVAL_GRAD, Q, dim),
    (qdata, :in, EVAL_NONE, Q, dim*(dim+1)÷2),
    (dv, :out, EVAL_GRAD, Q, dim),</code></pre><p>Note that the input array <code>du</code> has <a href="Globals.html#LibCEED.EvalMode"><code>EvalMode</code></a> <code>EVAL_GRAD</code>, meaning that this array stores the gradient of the trial function at each quadrature point. Therefore, at each quadrature point, <code>du</code> stores a vector of length <code>dim</code>, and so the shape of <code>du</code> is <code>(Q, dim)</code>. Similarly, the action of <span>$D$</span> is given by <span>$w \det(J) J^{-\intercal} J^{-1} \nabla u$</span>, which is also a vector of length <code>dim</code> at each quadrature point. This means that the output array <code>dv</code> also has shape <code>(Q, dim)</code>.</p><p>The geometric factors stored in <code>qdata</code> represent the symmetric matrix <span>$w \det(J) J^{-\intercal} J^{-1}$</span> evaluated at every quadrature point. In order to reduce data usage, instead of storing this data as a <span>$d \times d$</span> matrix, we use the fact that we know it is symmetric to only store <span>$d(d+1)/2$</span> entries, and the remaining entries we infer by symmetry. These entries are stored using the <a href="https://en.wikipedia.org/wiki/Voigt_notation">Voigt convention</a>. LibCEED.jl provides some <a href="Misc.html#LibCEED.getvoigt">utilities</a> for storing and extracting symmetric matrices stored in this fashion.</p><p>After the field specifications, we have the body of the Q-function:</p><pre><code class="language-julia">@inbounds @simd for i=1:Q
    dXdxdXdxT = getvoigt(@view(qdata[i,:]), CeedDim(dim))
    dui = SVector{dim}(@view(du[i,:]))
    dv[i,:] .= dXdxdXdxT*dui
end</code></pre><p>First, the matrix <span>$w \det(J) J^{-\intercal} J^{-1}$</span> is stored in the variable <code>dXdxdXdxT</code>. The symmetric entries of this matrix are accesed using <code>@view(qdata[i,:])</code>, which avoids allocations. <a href="Misc.html#LibCEED.getvoigt"><code>getvoigt</code></a> is used to convert from Voigt notation to a symmetric matrix, which returns a statically sized <code>SMatrix</code>. The version for the correct spatial dimension is selected using <code>CeedDim(dim)</code>, which allows for compile-time dispatch, since <code>dim</code> is a constant whose value is known as a constant when the Q-function is JIT compiled.</p><p>Then, the gradient of <span>$u$</span> at the given quadrature point is loaded as a fixed-size <code>SVector</code>. The result is placed into the output array, where the <a href="https://github.com/JuliaArrays/StaticArrays.jl"><code>StaticArrays.jl</code></a> package evaluates <code>dXdxdXdxT*dui</code> using an optimized matrix-vector product for small matrices (since their sizes are known statically).</p><h2 id="GPU-Kernels"><a class="docs-heading-anchor" href="#GPU-Kernels">GPU Kernels</a><a id="GPU-Kernels-1"></a><a class="docs-heading-anchor-permalink" href="#GPU-Kernels" title="Permalink"></a></h2><p>If the <code>Ceed</code> resource uses a CUDA backend, then the user Q-functions defined using <a href="QFunction.html#LibCEED.@interior_qf"><code>@interior_qf</code></a> are automatically compiled as CUDA kernels using <a href="https://github.com/JuliaGPU/CUDA.jl"><code>CUDA.jl</code></a>. Some Julia features are not available in GPU code (for example, dynamic dispatch), so if the Q-function is intended to be run on the GPU, the user should take care when defining the body of the user Q-function.</p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="C.html">« Low-level C interface</a><a class="docs-footer-nextpage" href="Examples.html">Examples »</a><div class="flexbox-break"></div><p class="footer-message">Powered by <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> and the <a href="https://julialang.org/">Julia Programming Language</a>.</p></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> on <span class="colophon-date" title="Wednesday 1 November 2023 06:18">Wednesday 1 November 2023</span>. Using Julia version 1.9.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>