1# Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2# the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3# reserved. See files LICENSE and NOTICE for details. 4# 5# This file is part of CEED, a collection of benchmarks, miniapps, software 6# libraries and APIs for efficient high-order finite element and spectral 7# element discretizations for exascale applications. For more information and 8# source code availability see http://github.com/ceed. 9# 10# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11# a collaborative effort of two U.S. Department of Energy organizations (Office 12# of Science and the National Nuclear Security Administration) responsible for 13# the planning and preparation of a capable exascale ecosystem, including 14# software, applications, hardware, advanced system engineering and early 15# testbed platforms, in support of the nation's exascale computing imperative. 16 17from _ceed_cffi import ffi, lib 18import sys 19import os 20import io 21import numpy as np 22import tempfile 23from abc import ABC 24from .ceed_vector import Vector 25from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1 26from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction 27from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction 28from .ceed_qfunctioncontext import QFunctionContext 29from .ceed_operator import Operator, CompositeOperator 30from .ceed_constants import * 31 32# ------------------------------------------------------------------------------ 33 34 35class Ceed(): 36 """Ceed: core components.""" 37 38 # Constructor 39 def __init__(self, resource="/cpu/self", on_error="store"): 40 # libCEED object 41 self._pointer = ffi.new("Ceed *") 42 43 # libCEED call 44 resourceAscii = ffi.new("char[]", resource.encode("ascii")) 45 os.environ["CEED_ERROR_HANDLER"] = "return" 46 err_code = lib.CeedInit(resourceAscii, self._pointer) 47 if err_code: 48 raise Exception("Error initializing backend resource: " + resource) 49 error_handlers = dict( 50 store="CeedErrorStore", 51 abort="CeedErrorAbort", 52 ) 53 lib.CeedSetErrorHandler( 54 self._pointer[0], ffi.addressof( 55 lib, error_handlers[on_error])) 56 57 # Representation 58 def __repr__(self): 59 return "<Ceed instance at " + hex(id(self)) + ">" 60 61 # String conversion for print() to stdout 62 def __str__(self): 63 """View a Ceed via print().""" 64 65 # libCEED call 66 with tempfile.NamedTemporaryFile() as key_file: 67 with open(key_file.name, 'r+') as stream_file: 68 stream = ffi.cast("FILE *", stream_file) 69 70 err_code = lib.CeedView(self._pointer[0], stream) 71 self._check_error(err_code) 72 73 stream_file.seek(0) 74 out_string = stream_file.read() 75 76 return out_string 77 78 # Error handler 79 def _check_error(self, err_code): 80 """Check return code and retrieve error message for non-zero code""" 81 if (err_code != lib.CEED_ERROR_SUCCESS): 82 message = ffi.new("char **") 83 lib.CeedGetErrorMessage(self._pointer[0], message) 84 raise Exception(ffi.string(message[0]).decode("UTF-8")) 85 86 # Get Resource 87 def get_resource(self): 88 """Get the full resource name for a Ceed context. 89 90 Returns: 91 resource: resource name""" 92 93 # libCEED call 94 resource = ffi.new("char **") 95 err_code = lib.CeedGetResource(self._pointer[0], resource) 96 self._check_error(err_code) 97 98 return ffi.string(resource[0]).decode("UTF-8") 99 100 # Preferred MemType 101 def get_preferred_memtype(self): 102 """Return Ceed preferred memory type. 103 104 Returns: 105 memtype: Ceed preferred memory type""" 106 107 # libCEED call 108 memtype = ffi.new("CeedMemType *", MEM_HOST) 109 err_code = lib.CeedGetPreferredMemType(self._pointer[0], memtype) 110 self._check_error(err_code) 111 112 return memtype[0] 113 114 # Convenience function to get CeedScalar type 115 def scalar_type(self): 116 """Return dtype string for CeedScalar. 117 118 Returns: 119 np_dtype: String naming numpy data type corresponding to CeedScalar""" 120 121 return scalar_types[lib.CEED_SCALAR_TYPE] 122 123 # --- Basis utility functions --- 124 125 # Gauss quadrature 126 def gauss_quadrature(self, q): 127 """Construct a Gauss-Legendre quadrature. 128 129 Args: 130 Q: number of quadrature points (integrates polynomials of 131 degree 2*Q-1 exactly) 132 133 Returns: 134 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 135 and array of length Q to hold the weights""" 136 137 # Setup arguments 138 qref1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 139 qweight1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 140 141 qref1d_pointer = ffi.new("CeedScalar *") 142 qref1d_pointer = ffi.cast( 143 "CeedScalar *", 144 qref1d.__array_interface__['data'][0]) 145 146 qweight1d_pointer = ffi.new("CeedScalar *") 147 qweight1d_pointer = ffi.cast( 148 "CeedScalar *", 149 qweight1d.__array_interface__['data'][0]) 150 151 # libCEED call 152 err_code = lib.CeedGaussQuadrature(q, qref1d_pointer, qweight1d_pointer) 153 self._check_error(err_code) 154 155 return qref1d, qweight1d 156 157 # Lobatto quadrature 158 def lobatto_quadrature(self, q): 159 """Construct a Gauss-Legendre-Lobatto quadrature. 160 161 Args: 162 q: number of quadrature points (integrates polynomials of 163 degree 2*Q-3 exactly) 164 165 Returns: 166 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 167 and array of length Q to hold the weights""" 168 169 # Setup arguments 170 qref1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 171 qref1d_pointer = ffi.new("CeedScalar *") 172 qref1d_pointer = ffi.cast( 173 "CeedScalar *", 174 qref1d.__array_interface__['data'][0]) 175 176 qweight1d = np.empty(q, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 177 qweight1d_pointer = ffi.new("CeedScalar *") 178 qweight1d_pointer = ffi.cast( 179 "CeedScalar *", 180 qweight1d.__array_interface__['data'][0]) 181 182 # libCEED call 183 err_code = lib.CeedLobattoQuadrature( 184 q, qref1d_pointer, qweight1d_pointer) 185 self._check_error(err_code) 186 187 return qref1d, qweight1d 188 189 # QR factorization 190 def qr_factorization(self, mat, tau, m, n): 191 """Return QR Factorization of a matrix. 192 193 Args: 194 ceed: Ceed context currently in use 195 *mat: Numpy array holding the row-major matrix to be factorized in place 196 *tau: Numpy array to hold the vector of lengt m of scaling factors 197 m: number of rows 198 n: numbef of columns""" 199 200 # Setup arguments 201 mat_pointer = ffi.new("CeedScalar *") 202 mat_pointer = ffi.cast( 203 "CeedScalar *", 204 mat.__array_interface__['data'][0]) 205 206 tau_pointer = ffi.new("CeedScalar *") 207 tau_pointer = ffi.cast( 208 "CeedScalar *", 209 tau.__array_interface__['data'][0]) 210 211 # libCEED call 212 err_code = lib.CeedQRFactorization( 213 self._pointer[0], mat_pointer, tau_pointer, m, n) 214 self._check_error(err_code) 215 216 return mat, tau 217 218 # Symmetric Schur decomposition 219 def symmetric_schur_decomposition(self, mat, n): 220 """Return symmetric Schur decomposition of a symmetric matrix 221 via symmetric QR factorization. 222 223 Args: 224 ceed: Ceed context currently in use 225 *mat: Numpy array holding the row-major matrix to be factorized in place 226 n: number of rows/columns 227 228 Returns: 229 lbda: Numpy array of length n holding eigenvalues""" 230 231 # Setup arguments 232 mat_pointer = ffi.new("CeedScalar *") 233 mat_pointer = ffi.cast( 234 "CeedScalar *", 235 mat.__array_interface__['data'][0]) 236 237 lbda = np.empty(n, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 238 l_pointer = ffi.new("CeedScalar *") 239 l_pointer = ffi.cast( 240 "CeedScalar *", 241 lbda.__array_interface__['data'][0]) 242 243 # libCEED call 244 err_code = lib.CeedSymmetricSchurDecomposition( 245 self._pointer[0], mat_pointer, l_pointer, n) 246 self._check_error(err_code) 247 248 return lbda 249 250 # Simultaneous Diagonalization 251 def simultaneous_diagonalization(self, matA, matB, n): 252 """Return Simultaneous Diagonalization of two matrices. 253 254 Args: 255 ceed: Ceed context currently in use 256 *matA: Numpy array holding the row-major matrix to be factorized with 257 eigenvalues 258 *matB: Numpy array holding the row-major matrix to be factorized to identity 259 n: number of rows/columns 260 261 Returns: 262 (x, lbda): Numpy array holding the row-major orthogonal matrix and 263 Numpy array holding the vector of length n of generalized 264 eigenvalues""" 265 266 # Setup arguments 267 matA_pointer = ffi.new("CeedScalar *") 268 matA_pointer = ffi.cast( 269 "CeedScalar *", 270 matA.__array_interface__['data'][0]) 271 272 matB_pointer = ffi.new("CeedScalar *") 273 matB_pointer = ffi.cast( 274 "CeedScalar *", 275 matB.__array_interface__['data'][0]) 276 277 lbda = np.empty(n, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 278 l_pointer = ffi.new("CeedScalar *") 279 l_pointer = ffi.cast( 280 "CeedScalar *", 281 lbda.__array_interface__['data'][0]) 282 283 x = np.empty(n * n, dtype=scalar_types[lib.CEED_SCALAR_TYPE]) 284 x_pointer = ffi.new("CeedScalar *") 285 x_pointer = ffi.cast("CeedScalar *", x.__array_interface__['data'][0]) 286 287 # libCEED call 288 err_code = lib.CeedSimultaneousDiagonalization(self._pointer[0], matA_pointer, matB_pointer, 289 x_pointer, l_pointer, n) 290 self._check_error(err_code) 291 292 return x, lbda 293 294 # --- libCEED objects --- 295 296 # CeedVector 297 def Vector(self, size): 298 """Ceed Vector: storing and manipulating vectors. 299 300 Args: 301 size: length of vector 302 303 Returns: 304 vector: Ceed Vector""" 305 306 return Vector(self, size) 307 308 # CeedElemRestriction 309 def ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, offsets, 310 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES): 311 """Ceed ElemRestriction: restriction from local vectors to elements. 312 313 Args: 314 nelem: number of elements described by the restriction 315 elemsize: size (number of nodes) per element 316 ncomp: number of field components per interpolation node 317 (1 for scalar fields) 318 compstride: Stride between components for the same L-vector "node". 319 Data for node i, component k can be found in the 320 L-vector at index [offsets[i] + k*compstride]. 321 lsize: The size of the L-vector. This vector may be larger than 322 the elements and fields given by this restriction. 323 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 324 ordered list of the offsets (into the input Ceed Vector) 325 for the unknowns corresponding to element i, where 326 0 <= i < nelem. All offsets must be in the range 327 [0, lsize - 1]. 328 **memtype: memory type of the offsets array, default CEED_MEM_HOST 329 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 330 331 Returns: 332 elemrestriction: Ceed ElemRestiction""" 333 334 return ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, 335 offsets, memtype=memtype, cmode=cmode) 336 337 def StridedElemRestriction(self, nelem, elemsize, ncomp, lsize, strides): 338 """Ceed Identity ElemRestriction: strided restriction from local vectors 339 to elements. 340 341 Args: 342 nelem: number of elements described by the restriction 343 elemsize: size (number of nodes) per element 344 ncomp: number of field components per interpolation node 345 (1 for scalar fields) 346 lsize: The size of the L-vector. This vector may be larger than 347 the elements and fields given by this restriction. 348 *strides: Array for strides between [nodes, components, elements]. 349 The data for node i, component j, element k in the 350 L-vector is given by 351 i*strides[0] + j*strides[1] + k*strides[2] 352 353 Returns: 354 elemrestriction: Ceed Strided ElemRestiction""" 355 356 return StridedElemRestriction( 357 self, nelem, elemsize, ncomp, lsize, strides) 358 359 def BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, compstride, 360 lsize, offsets, memtype=lib.CEED_MEM_HOST, 361 cmode=lib.CEED_COPY_VALUES): 362 """Ceed Blocked ElemRestriction: blocked restriction from local vectors 363 to elements. 364 365 Args: 366 nelem: number of elements described by the restriction 367 elemsize: size (number of nodes) per element 368 blksize: number of elements in a block 369 ncomp: number of field components per interpolation node 370 (1 for scalar fields) 371 lsize: The size of the L-vector. This vector may be larger than 372 the elements and fields given by this restriction. 373 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 374 ordered list of the offsets (into the input Ceed Vector) 375 for the unknowns corresponding to element i, where 376 0 <= i < nelem. All offsets must be in the range 377 [0, lsize - 1]. The backend will permute and pad this 378 array to the desired ordering for the blocksize, which is 379 typically given by the backend. The default reordering is 380 to interlace elements. 381 **memtype: memory type of the offsets array, default CEED_MEM_HOST 382 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 383 384 Returns: 385 elemrestriction: Ceed Blocked ElemRestiction""" 386 387 return BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, 388 compstride, lsize, offsets, 389 memtype=memtype, cmode=cmode) 390 391 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, ncomp, 392 lsize, strides): 393 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction 394 from local vectors to elements. 395 396 Args: 397 nelem: number of elements described in the restriction 398 elemsize: size (number of nodes) per element 399 blksize: number of elements in a block 400 ncomp: number of field components per interpolation node 401 (1 for scalar fields) 402 lsize: The size of the L-vector. This vector may be larger than 403 the elements and fields given by this restriction. 404 *strides: Array for strides between [nodes, components, elements]. 405 The data for node i, component j, element k in the 406 L-vector is given by 407 i*strides[0] + j*strides[1] + k*strides[2] 408 409 Returns: 410 elemrestriction: Ceed Strided ElemRestiction""" 411 412 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 413 ncomp, lsize, strides) 414 415 # CeedBasis 416 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 417 qref1d, qweight1d): 418 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 419 H^1 discretizations. 420 421 Args: 422 dim: topological dimension 423 ncomp: number of field components (1 for scalar fields) 424 P1d: number of nodes in one dimension 425 Q1d: number of quadrature points in one dimension 426 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix 427 expressing the values of nodal basis functions at 428 quadrature points 429 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix 430 expressing the derivatives of nodal basis functions at 431 quadrature points 432 *qref1d: Numpy array of length Q1d holding the locations of 433 quadrature points on the 1D reference element [-1, 1] 434 *qweight1d: Numpy array of length Q1d holding the quadrature 435 weights on the reference element 436 437 Returns: 438 basis: Ceed Basis""" 439 440 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 441 qref1d, qweight1d) 442 443 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 444 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange 445 basis objects for H^1 discretizations. 446 447 Args: 448 dim: topological dimension 449 ncomp: number of field components (1 for scalar fields) 450 P: number of Gauss-Lobatto nodes in one dimension. The 451 polynomial degree of the resulting Q_k element is k=P-1. 452 Q: number of quadrature points in one dimension 453 qmode: distribution of the Q quadrature points (affects order of 454 accuracy for the quadrature) 455 456 Returns: 457 basis: Ceed Basis""" 458 459 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 460 461 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 462 """Ceed H1 Basis: finite element non tensor-product basis for H^1 463 discretizations. 464 465 Args: 466 topo: topology of the element, e.g. hypercube, simplex, etc 467 ncomp: number of field components (1 for scalar fields) 468 nnodes: total number of nodes 469 nqpts: total number of quadrature points 470 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix 471 expressing the values of nodal basis functions at 472 quadrature points 473 *grad: Numpy array holding the row-major (nqpts * dim * nnodes) 474 matrix expressing the derivatives of nodal basis functions 475 at quadrature points 476 *qref: Numpy array of length (nqpts * dim) holding the locations of 477 quadrature points on the reference element [-1, 1] 478 *qweight: Numpy array of length nnodes holding the quadrature 479 weights on the reference element 480 481 Returns: 482 basis: Ceed Basis""" 483 484 return BasisH1(self, topo, ncomp, nnodes, nqpts, 485 interp, grad, qref, qweight) 486 487 # CeedQFunction 488 def QFunction(self, vlength, f, source): 489 """Ceed QFunction: point-wise operation at quadrature points for 490 evaluating volumetric terms. 491 492 Args: 493 vlength: vector length. Caller must ensure that number of quadrature 494 points is a multiple of vlength 495 f: ctypes function pointer to evaluate action at quadrature points 496 source: absolute path to source of QFunction, 497 "\\abs_path\\file.h:function_name 498 499 Returns: 500 qfunction: Ceed QFunction""" 501 502 return QFunction(self, vlength, f, source) 503 504 def QFunctionByName(self, name): 505 """Ceed QFunction By Name: point-wise operation at quadrature points 506 from a given gallery, for evaluating volumetric terms. 507 508 Args: 509 name: name of QFunction to use from gallery 510 511 Returns: 512 qfunction: Ceed QFunction By Name""" 513 514 return QFunctionByName(self, name) 515 516 def IdentityQFunction(self, size, inmode, outmode): 517 """Ceed Idenity QFunction: identity qfunction operation. 518 519 Args: 520 size: size of the qfunction fields 521 **inmode: CeedEvalMode for input to Ceed QFunction 522 **outmode: CeedEvalMode for output to Ceed QFunction 523 524 Returns: 525 qfunction: Ceed Identity QFunction""" 526 527 return IdentityQFunction(self, size, inmode, outmode) 528 529 def QFunctionContext(self): 530 """Ceed QFunction Context: stores Ceed QFunction user context data. 531 532 Returns: 533 userContext: Ceed QFunction Context""" 534 535 return QFunctionContext(self) 536 537 # CeedOperator 538 def Operator(self, qf, dqf=None, qdfT=None): 539 """Ceed Operator: composed FE-type operations on vectors. 540 541 Args: 542 qf: QFunction defining the action of the operator at quadrature 543 points 544 **dqf: QFunction defining the action of the Jacobian, default None 545 **dqfT: QFunction defining the action of the transpose of the 546 Jacobian, default None 547 548 Returns: 549 operator: Ceed Operator""" 550 551 return Operator(self, qf, dqf, qdfT) 552 553 def CompositeOperator(self): 554 """Composite Ceed Operator: composition of multiple CeedOperators. 555 556 Returns: 557 operator: Ceed Composite Operator""" 558 559 return CompositeOperator(self) 560 561 # Destructor 562 def __del__(self): 563 # libCEED call 564 lib.CeedDestroy(self._pointer) 565 566# ------------------------------------------------------------------------------ 567