1{ 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "metadata": {}, 6 "source": [ 7 "# libCEED for Python examples\n", 8 "\n", 9 "This is a tutorial to illustrate the main feautures of the Python interface for [libCEED](https://github.com/CEED/libCEED/), the low-level API library for efficient high-order discretization methods developed by the co-design [Center for Efficient Exascale Discretizations](https://ceed.exascaleproject.org/) (CEED) of the [Exascale Computing Project](https://www.exascaleproject.org/) (ECP).\n", 10 "\n", 11 "While libCEED's focus is on high-order finite/spectral element method implementations, the approach is mostly algebraic and thus applicable to other discretizations in factored form, as explained in the [user manual](https://libceed.org/)." 12 ] 13 }, 14 { 15 "cell_type": "markdown", 16 "metadata": {}, 17 "source": [ 18 "## Setting up libCEED for Python\n", 19 "\n", 20 "Install libCEED for Python by running" 21 ] 22 }, 23 { 24 "cell_type": "code", 25 "execution_count": null, 26 "metadata": {}, 27 "outputs": [], 28 "source": [ 29 "! python -m pip install libceed" 30 ] 31 }, 32 { 33 "cell_type": "markdown", 34 "metadata": {}, 35 "source": [ 36 "## CeedOperator\n", 37 "\n", 38 "Here we show some basic examples to illustrate the `libceed.Operator` class. In libCEED, a `libceed.Operator` defines the finite/spectral element operator associated to a `libceed.QFunction` (see [the API documentation](https://libceed.org/en/latest/libCEEDapi.html#finite-element-operator-decomposition))." 39 ] 40 }, 41 { 42 "cell_type": "markdown", 43 "metadata": {}, 44 "source": [ 45 "* In the following example, we create and apply a CeedOperator for the mass matrix in 1D. By applying this operator to a vector of 1's, we compute the length of this 1D domain, similar to Ex1-Volume in the [tutorial-6-shell tutorial](./tutorial-6-shell.ipynb)" 46 ] 47 }, 48 { 49 "cell_type": "code", 50 "execution_count": null, 51 "metadata": {}, 52 "outputs": [], 53 "source": [ 54 "import libceed\n", 55 "import numpy as np\n", 56 "\n", 57 "ceed = libceed.Ceed()\n", 58 "\n", 59 "num_elem = 15\n", 60 "p = 5\n", 61 "q = 8\n", 62 "num_x = num_elem + 1\n", 63 "num_u = num_elem*(p-1) + 1\n", 64 "\n", 65 "# Vectors\n", 66 "x = ceed.Vector(num_x)\n", 67 "x_array = np.zeros(num_x)\n", 68 "for i in range(num_x):\n", 69 " x_array[i] = i / (num_x - 1.0)\n", 70 "x.set_array(x_array, cmode=libceed.USE_POINTER)\n", 71 "\n", 72 "q_data = ceed.Vector(num_elem*q)\n", 73 "u = ceed.Vector(num_u)\n", 74 "v = ceed.Vector(num_u)\n", 75 "\n", 76 "# Restrictions\n", 77 "indices_x = np.zeros(num_x*2, dtype=\"int32\")\n", 78 "for i in range(num_x):\n", 79 " indices_x[2*i+0] = i\n", 80 " indices_x[2*i+1] = i+1\n", 81 "restriction_x = ceed.ElemRestriction(num_elem, 2, 1, 1, num_x, indices_x, cmode=libceed.USE_POINTER)\n", 82 "\n", 83 "indices_u = np.zeros(num_elem*p, dtype=\"int32\")\n", 84 "for i in range(num_elem):\n", 85 " for j in range(p):\n", 86 " indices_u[p*i+j] = i*(p-1) + j\n", 87 "restriction_u = ceed.ElemRestriction(num_elem, p, 1, 1, num_u, indices_u, cmode=libceed.USE_POINTER)\n", 88 "strides = np.array([1, q, q], dtype=\"int32\")\n", 89 "restriction_q_data = ceed.StridedElemRestriction(num_elem, q, 1, q*num_elem, strides)\n", 90 "\n", 91 "# Bases\n", 92 "basis_x = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS)\n", 93 "basis_u = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS)\n", 94 "\n", 95 "# QFunctions\n", 96 "qf_setup = ceed.QFunctionByName(\"Mass1DBuild\")\n", 97 "qf_mass = ceed.QFunctionByName(\"MassApply\")\n", 98 "\n", 99 "# Setup operator\n", 100 "op_setup = ceed.Operator(qf_setup)\n", 101 "op_setup.set_field(\"dx\", restriction_x, basis_x, libceed.VECTOR_ACTIVE)\n", 102 "op_setup.set_field(\"weights\", libceed.ELEMRESTRICTION_NONE, basis_x,\n", 103 " libceed.VECTOR_NONE)\n", 104 "op_setup.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE,\n", 105 " libceed.VECTOR_ACTIVE)\n", 106 "op_setup.check()\n", 107 "print('Setup operator: ', op_setup)\n", 108 "\n", 109 "# Mass operator\n", 110 "op_mass = ceed.Operator(qf_mass)\n", 111 "op_mass.set_field(\"u\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n", 112 "op_mass.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE, q_data)\n", 113 "op_mass.set_field(\"v\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n", 114 "op_mass.check()\n", 115 "print('Mass operator: ', op_mass)\n", 116 "\n", 117 "# Setup\n", 118 "op_setup.apply(x, q_data)\n", 119 "\n", 120 "# Apply mass matrix\n", 121 "u.set_value(1)\n", 122 "op_mass.apply(u, v)\n", 123 "\n", 124 "# Check\n", 125 "with v.array_read() as v_array:\n", 126 " print('The length of the domain is l = %4.2f'%np.sum(v_array))" 127 ] 128 }, 129 { 130 "cell_type": "markdown", 131 "metadata": {}, 132 "source": [ 133 "* In the next example, we create and apply a CeedOperator for the Poisson operator in 1D. By applying this operator to a vector with a linear function, we compute the 'surface area' of this 1D domain, similar to Ex2-Surface in the [tutorial-6-shell tutorial](./tutorial-6-shell.ipynb)" 134 ] 135 }, 136 { 137 "cell_type": "code", 138 "execution_count": null, 139 "metadata": {}, 140 "outputs": [], 141 "source": [ 142 "import libceed\n", 143 "import numpy as np\n", 144 "\n", 145 "ceed = libceed.Ceed()\n", 146 "\n", 147 "num_elem = 15\n", 148 "p = 5\n", 149 "q = 8\n", 150 "num_x = num_elem + 1\n", 151 "num_u = num_elem*(p-1) + 1\n", 152 "\n", 153 "# Vectors\n", 154 "x = ceed.Vector(num_x)\n", 155 "x_array = np.zeros(num_x)\n", 156 "for i in range(num_x):\n", 157 " x_array[i] = i / (num_x - 1.0)\n", 158 "x.set_array(x_array, cmode=libceed.USE_POINTER)\n", 159 "\n", 160 "q_data = ceed.Vector(num_elem*q)\n", 161 "u = ceed.Vector(num_u)\n", 162 "v = ceed.Vector(num_u)\n", 163 "\n", 164 "# Restrictions\n", 165 "indices_x = np.zeros(num_x*2, dtype=\"int32\")\n", 166 "for i in range(num_x):\n", 167 " indices_x[2*i+0] = i\n", 168 " indices_x[2*i+1] = i+1\n", 169 "restriction_x = ceed.ElemRestriction(num_elem, 2, 1, 1, num_x, indices_x, cmode=libceed.USE_POINTER)\n", 170 "\n", 171 "indices_u = np.zeros(num_elem*p, dtype=\"int32\")\n", 172 "for i in range(num_elem):\n", 173 " for j in range(p):\n", 174 " indices_u[p*i+j] = i*(p-1) + j\n", 175 "restriction_u = ceed.ElemRestriction(num_elem, p, 1, 1, num_u, indices_u, cmode=libceed.USE_POINTER)\n", 176 "strides = np.array([1, q, q], dtype=\"int32\")\n", 177 "restriction_q_data = ceed.StridedElemRestriction(num_elem, q, 1, q*num_elem, strides)\n", 178 "\n", 179 "# Bases\n", 180 "basis_x = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS)\n", 181 "basis_u = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS)\n", 182 "\n", 183 "# QFunctions\n", 184 "qf_setup = ceed.QFunctionByName(\"Poisson1DBuild\")\n", 185 "qf_mass = ceed.QFunctionByName(\"Poisson1DApply\")\n", 186 "\n", 187 "# Setup operator\n", 188 "op_setup = ceed.Operator(qf_setup)\n", 189 "op_setup.set_field(\"dx\", restriction_x, basis_x, libceed.VECTOR_ACTIVE)\n", 190 "op_setup.set_field(\"weights\", libceed.ELEMRESTRICTION_NONE, basis_x,\n", 191 " libceed.VECTOR_NONE)\n", 192 "op_setup.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE,\n", 193 " libceed.VECTOR_ACTIVE)\n", 194 "op_setup.check()\n", 195 "print('Setup operator: ', op_setup)\n", 196 "\n", 197 "# Poisson operator\n", 198 "op_poisson = ceed.Operator(qf_mass)\n", 199 "op_poisson.set_field(\"du\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n", 200 "op_poisson.set_field(\"qdata\", restriction_q_data, libceed.BASIS_NONE, q_data)\n", 201 "op_poisson.set_field(\"dv\", restriction_u, basis_u, libceed.VECTOR_ACTIVE)\n", 202 "op_poisson.check()\n", 203 "print('Poisson operator: ', op_poisson)\n", 204 "\n", 205 "# Setup\n", 206 "op_setup.apply(x, q_data)\n", 207 "\n", 208 "# Apply Poisson operator\n", 209 "with u.array_write() as u_array:\n", 210 " [points, _] = ceed.lobatto_quadrature(p)\n", 211 " for elem in range(num_elem):\n", 212 " for point in range(p):\n", 213 " u_array[elem * (p - 1) + point] = (1.0 + 2.0 * elem + points[point])/(2.0 * num_elem)\n", 214 "op_poisson.apply(u, v)\n", 215 "\n", 216 "# Check\n", 217 "with v.array_read() as v_array:\n", 218 " print('The surface area of the domain is dl = %4.2f'%np.sum(abs(v_array)))" 219 ] 220 } 221 ], 222 "metadata": { 223 "kernelspec": { 224 "display_name": "Python 3 (ipykernel)", 225 "language": "python", 226 "name": "python3" 227 }, 228 "language_info": { 229 "codemirror_mode": { 230 "name": "ipython", 231 "version": 3 232 }, 233 "file_extension": ".py", 234 "mimetype": "text/x-python", 235 "name": "python", 236 "nbconvert_exporter": "python", 237 "pygments_lexer": "ipython3", 238 "version": "3.13.2" 239 } 240 }, 241 "nbformat": 4, 242 "nbformat_minor": 4 243} 244