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.readthedocs.io/)." 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 "## CeedBasis\n", 37 "\n", 38 "Here we show some basic examples to illustrate the `libceed.Basis` class. In libCEED, a `libceed.Basis` defines the finite element basis and associated quadrature rule (see [the API documentation](https://libceed.readthedocs.io/en/latest/libCEEDapi.html#finite-element-operator-decomposition))." 39 ] 40 }, 41 { 42 "cell_type": "markdown", 43 "metadata": {}, 44 "source": [ 45 "First we declare some auxiliary functions needed in the following examples" 46 ] 47 }, 48 { 49 "cell_type": "code", 50 "execution_count": null, 51 "metadata": {}, 52 "outputs": [], 53 "source": [ 54 "%matplotlib inline\n", 55 "import numpy as np\n", 56 "import matplotlib.pyplot as plt\n", 57 "plt.style.use('ggplot')\n", 58 "\n", 59 "def eval(dim, x):\n", 60 " result, center = 1, 0.1\n", 61 " for d in range(dim):\n", 62 " result *= np.tanh(x[d] - center)\n", 63 " center += 0.1\n", 64 " return result\n", 65 "\n", 66 "def feval(x1, x2):\n", 67 " return x1*x1 + x2*x2 + x1*x2 + 1\n", 68 "\n", 69 "def dfeval(x1, x2):\n", 70 " return 2*x1 + x2" 71 ] 72 }, 73 { 74 "cell_type": "markdown", 75 "metadata": {}, 76 "source": [ 77 "## $H^1$ Lagrange bases in 1D\n", 78 "\n", 79 "The Lagrange interpolation nodes are at the Gauss-Lobatto points, so interpolation to Gauss-Lobatto quadrature points is the identity." 80 ] 81 }, 82 { 83 "cell_type": "code", 84 "execution_count": null, 85 "metadata": {}, 86 "outputs": [], 87 "source": [ 88 "import libceed\n", 89 "\n", 90 "ceed = libceed.Ceed()\n", 91 "\n", 92 "b = ceed.BasisTensorH1Lagrange(\n", 93 " dim=1, # topological dimension\n", 94 " ncomp=1, # number of components\n", 95 " P=4, # number of basis functions (nodes) per dimension\n", 96 " Q=4, # number of quadrature points per dimension\n", 97 " qmode=libceed.GAUSS_LOBATTO)\n", 98 "print(b)" 99 ] 100 }, 101 { 102 "cell_type": "markdown", 103 "metadata": {}, 104 "source": [ 105 "Although a `libceed.Basis` is fully discrete, we can use the Lagrange construction to extend the basis to continuous functions by applying `EVAL_INTERP` to the identity. This is the Vandermonde matrix of the continuous basis." 106 ] 107 }, 108 { 109 "cell_type": "code", 110 "execution_count": null, 111 "metadata": {}, 112 "outputs": [], 113 "source": [ 114 "P = b.get_num_nodes()\n", 115 "nviz = 50\n", 116 "bviz = ceed.BasisTensorH1Lagrange(1, 1, P, nviz, libceed.GAUSS_LOBATTO)\n", 117 "\n", 118 "# Construct P \"elements\" with one node activated\n", 119 "I = ceed.Vector(P * P)\n", 120 "with I.array(P, P) as x:\n", 121 " x[...] = np.eye(P)\n", 122 "\n", 123 "Bvander = ceed.Vector(P * nviz)\n", 124 "bviz.apply(4, libceed.EVAL_INTERP, I, Bvander)\n", 125 "\n", 126 "qviz, _weight = ceed.lobatto_quadrature(nviz)\n", 127 "with Bvander.array_read(nviz, P) as B:\n", 128 " plt.plot(qviz, B)\n", 129 "\n", 130 "# Mark tho Lobatto nodes\n", 131 "qb, _weight = ceed.lobatto_quadrature(P)\n", 132 "plt.plot(qb, 0*qb, 'ok');" 133 ] 134 }, 135 { 136 "cell_type": "markdown", 137 "metadata": {}, 138 "source": [ 139 "In contrast, the Gauss quadrature points are not collocated, and thus all basis functions are generally nonzero at every quadrature point." 140 ] 141 }, 142 { 143 "cell_type": "code", 144 "execution_count": null, 145 "metadata": {}, 146 "outputs": [], 147 "source": [ 148 "b = ceed.BasisTensorH1Lagrange(1, 1, 4, 4, libceed.GAUSS)\n", 149 "print(b)\n", 150 "\n", 151 "with Bvander.array_read(nviz, P) as B:\n", 152 " plt.plot(qviz, B)\n", 153 "# Mark tho Gauss quadrature points\n", 154 "qb, _weight = ceed.gauss_quadrature(P)\n", 155 "plt.plot(qb, 0*qb, 'ok');" 156 ] 157 }, 158 { 159 "cell_type": "markdown", 160 "metadata": {}, 161 "source": [ 162 "Although the underlying functions are not an intrinsic property of a `libceed.Basis` in libCEED, the sizes are.\n", 163 "Here, we create a 3D tensor product element with more quadrature points than Lagrange interpolation nodes." 164 ] 165 }, 166 { 167 "cell_type": "code", 168 "execution_count": null, 169 "metadata": {}, 170 "outputs": [], 171 "source": [ 172 "b = ceed.BasisTensorH1Lagrange(3, 1, 4, 5, libceed.GAUSS_LOBATTO)\n", 173 "\n", 174 "p = b.get_num_nodes()\n", 175 "print('p =', p)\n", 176 "\n", 177 "q = b.get_num_quadrature_points()\n", 178 "print('q =', q)" 179 ] 180 }, 181 { 182 "cell_type": "markdown", 183 "metadata": {}, 184 "source": [ 185 "* In the following example, we demonstrate the application of an interpolatory basis in multiple dimensions" 186 ] 187 }, 188 { 189 "cell_type": "code", 190 "execution_count": null, 191 "metadata": {}, 192 "outputs": [], 193 "source": [ 194 "for dim in range(1, 4):\n", 195 " Q = 4\n", 196 " Qdim = Q**dim\n", 197 " Xdim = 2**dim\n", 198 " x = np.empty(Xdim*dim, dtype=\"float64\")\n", 199 " uq = np.empty(Qdim, dtype=\"float64\")\n", 200 "\n", 201 " for d in range(dim):\n", 202 " for i in range(Xdim):\n", 203 " x[d*Xdim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n", 204 "\n", 205 " X = ceed.Vector(Xdim*dim)\n", 206 " X.set_array(x, cmode=libceed.USE_POINTER)\n", 207 " Xq = ceed.Vector(Qdim*dim)\n", 208 " Xq.set_value(0)\n", 209 " U = ceed.Vector(Qdim)\n", 210 " U.set_value(0)\n", 211 " Uq = ceed.Vector(Qdim)\n", 212 "\n", 213 " bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)\n", 214 " bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)\n", 215 "\n", 216 " bxl.apply(1, libceed.EVAL_INTERP, X, Xq)\n", 217 "\n", 218 " with Xq.array_read() as xq:\n", 219 " for i in range(Qdim):\n", 220 " xx = np.empty(dim, dtype=\"float64\")\n", 221 " for d in range(dim):\n", 222 " xx[d] = xq[d*Qdim + i]\n", 223 " uq[i] = eval(dim, xx)\n", 224 "\n", 225 " Uq.set_array(uq, cmode=libceed.USE_POINTER)\n", 226 "\n", 227 " # This operation is the identity because the quadrature is collocated\n", 228 " bul.T.apply(1, libceed.EVAL_INTERP, Uq, U)\n", 229 "\n", 230 " bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)\n", 231 " bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)\n", 232 "\n", 233 " bxg.apply(1, libceed.EVAL_INTERP, X, Xq)\n", 234 " bug.apply(1, libceed.EVAL_INTERP, U, Uq)\n", 235 "\n", 236 " with Xq.array_read() as xq, Uq.array_read() as u:\n", 237 " #print('xq =', xq)\n", 238 " #print('u =', u)\n", 239 " if dim == 2:\n", 240 " # Default ordering is contiguous in x direction, but\n", 241 " # pyplot expects meshgrid convention, which is transposed.\n", 242 " x, y = xq.reshape(2, Q, Q).transpose(0, 2, 1)\n", 243 " plt.scatter(x, y, c=np.array(u).reshape(Q, Q))\n", 244 " plt.xlim(-1, 1)\n", 245 " plt.ylim(-1, 1)\n", 246 " plt.colorbar(label='u')" 247 ] 248 }, 249 { 250 "cell_type": "markdown", 251 "metadata": {}, 252 "source": [ 253 "* In the following example, we demonstrate the application of the gradient of the shape functions in multiple dimensions" 254 ] 255 }, 256 { 257 "cell_type": "code", 258 "execution_count": null, 259 "metadata": {}, 260 "outputs": [], 261 "source": [ 262 "for dim in range (1, 4):\n", 263 " P, Q = 8, 10\n", 264 " Pdim = P**dim\n", 265 " Qdim = Q**dim\n", 266 " Xdim = 2**dim\n", 267 " sum1 = sum2 = 0\n", 268 " x = np.empty(Xdim*dim, dtype=\"float64\")\n", 269 " u = np.empty(Pdim, dtype=\"float64\")\n", 270 "\n", 271 " for d in range(dim):\n", 272 " for i in range(Xdim):\n", 273 " x[d*Xdim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1\n", 274 "\n", 275 " X = ceed.Vector(Xdim*dim)\n", 276 " X.set_array(x, cmode=libceed.USE_POINTER)\n", 277 " Xq = ceed.Vector(Pdim*dim)\n", 278 " Xq.set_value(0)\n", 279 " U = ceed.Vector(Pdim)\n", 280 " Uq = ceed.Vector(Qdim*dim)\n", 281 " Uq.set_value(0)\n", 282 " Ones = ceed.Vector(Qdim*dim)\n", 283 " Ones.set_value(1)\n", 284 " Gtposeones = ceed.Vector(Pdim)\n", 285 " Gtposeones.set_value(0)\n", 286 "\n", 287 " # Get function values at quadrature points\n", 288 " bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, P, libceed.GAUSS_LOBATTO)\n", 289 " bxl.apply(1, libceed.EVAL_INTERP, X, Xq)\n", 290 "\n", 291 " with Xq.array_read() as xq:\n", 292 " for i in range(Pdim):\n", 293 " xx = np.empty(dim, dtype=\"float64\")\n", 294 " for d in range(dim):\n", 295 " xx[d] = xq[d*Pdim + i]\n", 296 " u[i] = eval(dim, xx)\n", 297 "\n", 298 " U.set_array(u, cmode=libceed.USE_POINTER)\n", 299 "\n", 300 " # Calculate G u at quadrature points, G' * 1 at dofs\n", 301 " bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)\n", 302 " bug.apply(1, libceed.EVAL_GRAD, U, Uq)\n", 303 " bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones)\n", 304 "\n", 305 " # Check if 1' * G * u = u' * (G' * 1)\n", 306 " with Gtposeones.array_read() as gtposeones, Uq.array_read() as uq:\n", 307 " for i in range(Pdim):\n", 308 " sum1 += gtposeones[i]*u[i]\n", 309 " for i in range(dim*Qdim):\n", 310 " sum2 += uq[i]\n", 311 "\n", 312 " # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n", 313 " print('1T * G * u - uT * (GT * 1) =', np.abs(sum1 - sum2))" 314 ] 315 }, 316 { 317 "cell_type": "markdown", 318 "metadata": {}, 319 "source": [ 320 "### Advanced topics" 321 ] 322 }, 323 { 324 "cell_type": "markdown", 325 "metadata": {}, 326 "source": [ 327 "* In the following example, we demonstrate the QR factorization of a basis matrix.\n", 328 "The representation is similar to LAPACK's [`dgeqrf`](https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html#ga3766ea903391b5cf9008132f7440ec7b), with elementary reflectors in the lower triangular block, scaled by `tau`." 329 ] 330 }, 331 { 332 "cell_type": "code", 333 "execution_count": null, 334 "metadata": {}, 335 "outputs": [], 336 "source": [ 337 "qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype=\"float64\")\n", 338 "tau = np.empty(3, dtype=\"float64\")\n", 339 "\n", 340 "qr, tau = ceed.qr_factorization(qr, tau, 4, 3)\n", 341 "\n", 342 "print('qr =')\n", 343 "print(qr.reshape(4, 3))\n", 344 "\n", 345 "print('tau =')\n", 346 "print(tau)" 347 ] 348 }, 349 { 350 "cell_type": "markdown", 351 "metadata": {}, 352 "source": [ 353 "* In the following example, we demonstrate the symmetric Schur decomposition of a basis matrix" 354 ] 355 }, 356 { 357 "cell_type": "code", 358 "execution_count": null, 359 "metadata": {}, 360 "outputs": [], 361 "source": [ 362 "A = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n", 363 " 0.0745459, 1., 0.16666509, -0.07448852,\n", 364 " -0.07448852, 0.16666509, 1., 0.0745459,\n", 365 " 0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n", 366 "\n", 367 "lam = ceed.symmetric_schur_decomposition(A, 4)\n", 368 "\n", 369 "print(\"Q =\")\n", 370 "for i in range(4):\n", 371 " for j in range(4):\n", 372 " if A[j+4*i] <= 1E-14 and A[j+4*i] >= -1E-14:\n", 373 " A[j+4*i] = 0\n", 374 " print(\"%12.8f\"%A[j+4*i])\n", 375 "\n", 376 "print(\"lambda =\")\n", 377 "for i in range(4):\n", 378 " if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n", 379 " lam[i] = 0\n", 380 " print(\"%12.8f\"%lam[i])" 381 ] 382 }, 383 { 384 "cell_type": "markdown", 385 "metadata": {}, 386 "source": [ 387 "* In the following example, we demonstrate the simultaneous diagonalization of a basis matrix" 388 ] 389 }, 390 { 391 "cell_type": "code", 392 "execution_count": null, 393 "metadata": {}, 394 "outputs": [], 395 "source": [ 396 "M = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n", 397 " 0.0745459, 1., 0.16666509, -0.07448852,\n", 398 " -0.07448852, 0.16666509, 1., 0.0745459,\n", 399 " 0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n", 400 "K = np.array([3.03344425, -3.41501767, 0.49824435, -0.11667092,\n", 401 " -3.41501767, 5.83354662, -2.9167733, 0.49824435,\n", 402 " 0.49824435, -2.9167733, 5.83354662, -3.41501767,\n", 403 " -0.11667092, 0.49824435, -3.41501767, 3.03344425], dtype=\"float64\")\n", 404 "\n", 405 "x, lam = ceed.simultaneous_diagonalization(K, M, 4)\n", 406 "\n", 407 "print(\"x =\")\n", 408 "for i in range(4):\n", 409 " for j in range(4):\n", 410 " if x[j+4*i] <= 1E-14 and x[j+4*i] >= -1E-14:\n", 411 " x[j+4*i] = 0\n", 412 " print(\"%12.8f\"%x[j+4*i])\n", 413 "\n", 414 "print(\"lambda =\")\n", 415 "for i in range(4):\n", 416 " if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n", 417 " lam[i] = 0\n", 418 " print(\"%12.8f\"%lam[i])" 419 ] 420 } 421 ], 422 "metadata": { 423 "kernelspec": { 424 "display_name": "Python 3", 425 "language": "python", 426 "name": "python3" 427 }, 428 "language_info": { 429 "codemirror_mode": { 430 "name": "ipython", 431 "version": 3 432 }, 433 "file_extension": ".py", 434 "mimetype": "text/x-python", 435 "name": "python", 436 "nbconvert_exporter": "python", 437 "pygments_lexer": "ipython3", 438 "version": "3.8.5" 439 } 440 }, 441 "nbformat": 4, 442 "nbformat_minor": 4 443} 444