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 # --- Basis utility functions --- 115 116 # Gauss quadrature 117 def gauss_quadrature(self, q): 118 """Construct a Gauss-Legendre quadrature. 119 120 Args: 121 Q: number of quadrature points (integrates polynomials of 122 degree 2*Q-1 exactly) 123 124 Returns: 125 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 126 and array of length Q to hold the weights""" 127 128 # Setup arguments 129 qref1d = np.empty(q, dtype="float64") 130 qweight1d = np.empty(q, dtype="float64") 131 132 qref1d_pointer = ffi.new("CeedScalar *") 133 qref1d_pointer = ffi.cast( 134 "CeedScalar *", 135 qref1d.__array_interface__['data'][0]) 136 137 qweight1d_pointer = ffi.new("CeedScalar *") 138 qweight1d_pointer = ffi.cast( 139 "CeedScalar *", 140 qweight1d.__array_interface__['data'][0]) 141 142 # libCEED call 143 err_code = lib.CeedGaussQuadrature(q, qref1d_pointer, qweight1d_pointer) 144 self._check_error(err_code) 145 146 return qref1d, qweight1d 147 148 # Lobatto quadrature 149 def lobatto_quadrature(self, q): 150 """Construct a Gauss-Legendre-Lobatto quadrature. 151 152 Args: 153 q: number of quadrature points (integrates polynomials of 154 degree 2*Q-3 exactly) 155 156 Returns: 157 (qref1d, qweight1d): array of length Q to hold the abscissa on [-1, 1] 158 and array of length Q to hold the weights""" 159 160 # Setup arguments 161 qref1d = np.empty(q, dtype="float64") 162 qref1d_pointer = ffi.new("CeedScalar *") 163 qref1d_pointer = ffi.cast( 164 "CeedScalar *", 165 qref1d.__array_interface__['data'][0]) 166 167 qweight1d = np.empty(q, dtype="float64") 168 qweight1d_pointer = ffi.new("CeedScalar *") 169 qweight1d_pointer = ffi.cast( 170 "CeedScalar *", 171 qweight1d.__array_interface__['data'][0]) 172 173 # libCEED call 174 err_code = lib.CeedLobattoQuadrature( 175 q, qref1d_pointer, qweight1d_pointer) 176 self._check_error(err_code) 177 178 return qref1d, qweight1d 179 180 # QR factorization 181 def qr_factorization(self, mat, tau, m, n): 182 """Return QR Factorization of a matrix. 183 184 Args: 185 ceed: Ceed context currently in use 186 *mat: Numpy array holding the row-major matrix to be factorized in place 187 *tau: Numpy array to hold the vector of lengt m of scaling factors 188 m: number of rows 189 n: numbef of columns""" 190 191 # Setup arguments 192 mat_pointer = ffi.new("CeedScalar *") 193 mat_pointer = ffi.cast( 194 "CeedScalar *", 195 mat.__array_interface__['data'][0]) 196 197 tau_pointer = ffi.new("CeedScalar *") 198 tau_pointer = ffi.cast( 199 "CeedScalar *", 200 tau.__array_interface__['data'][0]) 201 202 # libCEED call 203 err_code = lib.CeedQRFactorization( 204 self._pointer[0], mat_pointer, tau_pointer, m, n) 205 self._check_error(err_code) 206 207 return mat, tau 208 209 # Symmetric Schur decomposition 210 def symmetric_schur_decomposition(self, mat, n): 211 """Return symmetric Schur decomposition of a symmetric matrix 212 via symmetric QR factorization. 213 214 Args: 215 ceed: Ceed context currently in use 216 *mat: Numpy array holding the row-major matrix to be factorized in place 217 n: number of rows/columns 218 219 Returns: 220 lbda: Numpy array of length n holding eigenvalues""" 221 222 # Setup arguments 223 mat_pointer = ffi.new("CeedScalar *") 224 mat_pointer = ffi.cast( 225 "CeedScalar *", 226 mat.__array_interface__['data'][0]) 227 228 lbda = np.empty(n, dtype="float64") 229 l_pointer = ffi.new("CeedScalar *") 230 l_pointer = ffi.cast( 231 "CeedScalar *", 232 lbda.__array_interface__['data'][0]) 233 234 # libCEED call 235 err_code = lib.CeedSymmetricSchurDecomposition( 236 self._pointer[0], mat_pointer, l_pointer, n) 237 self._check_error(err_code) 238 239 return lbda 240 241 # Simultaneous Diagonalization 242 def simultaneous_diagonalization(self, matA, matB, n): 243 """Return Simultaneous Diagonalization of two matrices. 244 245 Args: 246 ceed: Ceed context currently in use 247 *matA: Numpy array holding the row-major matrix to be factorized with 248 eigenvalues 249 *matB: Numpy array holding the row-major matrix to be factorized to identity 250 n: number of rows/columns 251 252 Returns: 253 (x, lbda): Numpy array holding the row-major orthogonal matrix and 254 Numpy array holding the vector of length n of generalized 255 eigenvalues""" 256 257 # Setup arguments 258 matA_pointer = ffi.new("CeedScalar *") 259 matA_pointer = ffi.cast( 260 "CeedScalar *", 261 matA.__array_interface__['data'][0]) 262 263 matB_pointer = ffi.new("CeedScalar *") 264 matB_pointer = ffi.cast( 265 "CeedScalar *", 266 matB.__array_interface__['data'][0]) 267 268 lbda = np.empty(n, dtype="float64") 269 l_pointer = ffi.new("CeedScalar *") 270 l_pointer = ffi.cast( 271 "CeedScalar *", 272 lbda.__array_interface__['data'][0]) 273 274 x = np.empty(n * n, dtype="float64") 275 x_pointer = ffi.new("CeedScalar *") 276 x_pointer = ffi.cast("CeedScalar *", x.__array_interface__['data'][0]) 277 278 # libCEED call 279 err_code = lib.CeedSimultaneousDiagonalization(self._pointer[0], matA_pointer, matB_pointer, 280 x_pointer, l_pointer, n) 281 self._check_error(err_code) 282 283 return x, lbda 284 285 # --- libCEED objects --- 286 287 # CeedVector 288 def Vector(self, size): 289 """Ceed Vector: storing and manipulating vectors. 290 291 Args: 292 size: length of vector 293 294 Returns: 295 vector: Ceed Vector""" 296 297 return Vector(self, size) 298 299 # CeedElemRestriction 300 def ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, offsets, 301 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES): 302 """Ceed ElemRestriction: restriction from local vectors to elements. 303 304 Args: 305 nelem: number of elements described by the restriction 306 elemsize: size (number of nodes) per element 307 ncomp: number of field components per interpolation node 308 (1 for scalar fields) 309 compstride: Stride between components for the same L-vector "node". 310 Data for node i, component k can be found in the 311 L-vector at index [offsets[i] + k*compstride]. 312 lsize: The size of the L-vector. This vector may be larger than 313 the elements and fields given by this restriction. 314 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 315 ordered list of the offsets (into the input Ceed Vector) 316 for the unknowns corresponding to element i, where 317 0 <= i < nelem. All offsets must be in the range 318 [0, lsize - 1]. 319 **memtype: memory type of the offsets array, default CEED_MEM_HOST 320 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 321 322 Returns: 323 elemrestriction: Ceed ElemRestiction""" 324 325 return ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, 326 offsets, memtype=memtype, cmode=cmode) 327 328 def StridedElemRestriction(self, nelem, elemsize, ncomp, lsize, strides): 329 """Ceed Identity ElemRestriction: strided restriction from local vectors 330 to elements. 331 332 Args: 333 nelem: number of elements described by the restriction 334 elemsize: size (number of nodes) per element 335 ncomp: number of field components per interpolation node 336 (1 for scalar fields) 337 lsize: The size of the L-vector. This vector may be larger than 338 the elements and fields given by this restriction. 339 *strides: Array for strides between [nodes, components, elements]. 340 The data for node i, component j, element k in the 341 L-vector is given by 342 i*strides[0] + j*strides[1] + k*strides[2] 343 344 Returns: 345 elemrestriction: Ceed Strided ElemRestiction""" 346 347 return StridedElemRestriction( 348 self, nelem, elemsize, ncomp, lsize, strides) 349 350 def BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, compstride, 351 lsize, offsets, memtype=lib.CEED_MEM_HOST, 352 cmode=lib.CEED_COPY_VALUES): 353 """Ceed Blocked ElemRestriction: blocked restriction from local vectors 354 to elements. 355 356 Args: 357 nelem: number of elements described by the restriction 358 elemsize: size (number of nodes) per element 359 blksize: number of elements in a block 360 ncomp: number of field components per interpolation node 361 (1 for scalar fields) 362 lsize: The size of the L-vector. This vector may be larger than 363 the elements and fields given by this restriction. 364 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 365 ordered list of the offsets (into the input Ceed Vector) 366 for the unknowns corresponding to element i, where 367 0 <= i < nelem. All offsets must be in the range 368 [0, lsize - 1]. The backend will permute and pad this 369 array to the desired ordering for the blocksize, which is 370 typically given by the backend. The default reordering is 371 to interlace elements. 372 **memtype: memory type of the offsets array, default CEED_MEM_HOST 373 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 374 375 Returns: 376 elemrestriction: Ceed Blocked ElemRestiction""" 377 378 return BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, 379 compstride, lsize, offsets, 380 memtype=memtype, cmode=cmode) 381 382 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, ncomp, 383 lsize, strides): 384 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction 385 from local vectors to elements. 386 387 Args: 388 nelem: number of elements described in the restriction 389 elemsize: size (number of nodes) per element 390 blksize: number of elements in a block 391 ncomp: number of field components per interpolation node 392 (1 for scalar fields) 393 lsize: The size of the L-vector. This vector may be larger than 394 the elements and fields given by this restriction. 395 *strides: Array for strides between [nodes, components, elements]. 396 The data for node i, component j, element k in the 397 L-vector is given by 398 i*strides[0] + j*strides[1] + k*strides[2] 399 400 Returns: 401 elemrestriction: Ceed Strided ElemRestiction""" 402 403 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 404 ncomp, lsize, strides) 405 406 # CeedBasis 407 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 408 qref1d, qweight1d): 409 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 410 H^1 discretizations. 411 412 Args: 413 dim: topological dimension 414 ncomp: number of field components (1 for scalar fields) 415 P1d: number of nodes in one dimension 416 Q1d: number of quadrature points in one dimension 417 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix 418 expressing the values of nodal basis functions at 419 quadrature points 420 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix 421 expressing the derivatives of nodal basis functions at 422 quadrature points 423 *qref1d: Numpy array of length Q1d holding the locations of 424 quadrature points on the 1D reference element [-1, 1] 425 *qweight1d: Numpy array of length Q1d holding the quadrature 426 weights on the reference element 427 428 Returns: 429 basis: Ceed Basis""" 430 431 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 432 qref1d, qweight1d) 433 434 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 435 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange 436 basis objects for H^1 discretizations. 437 438 Args: 439 dim: topological dimension 440 ncomp: number of field components (1 for scalar fields) 441 P: number of Gauss-Lobatto nodes in one dimension. The 442 polynomial degree of the resulting Q_k element is k=P-1. 443 Q: number of quadrature points in one dimension 444 qmode: distribution of the Q quadrature points (affects order of 445 accuracy for the quadrature) 446 447 Returns: 448 basis: Ceed Basis""" 449 450 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 451 452 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 453 """Ceed H1 Basis: finite element non tensor-product basis for H^1 454 discretizations. 455 456 Args: 457 topo: topology of the element, e.g. hypercube, simplex, etc 458 ncomp: number of field components (1 for scalar fields) 459 nnodes: total number of nodes 460 nqpts: total number of quadrature points 461 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix 462 expressing the values of nodal basis functions at 463 quadrature points 464 *grad: Numpy array holding the row-major (nqpts * dim * nnodes) 465 matrix expressing the derivatives of nodal basis functions 466 at quadrature points 467 *qref: Numpy array of length (nqpts * dim) holding the locations of 468 quadrature points on the reference element [-1, 1] 469 *qweight: Numpy array of length nnodes holding the quadrature 470 weights on the reference element 471 472 Returns: 473 basis: Ceed Basis""" 474 475 return BasisH1(self, topo, ncomp, nnodes, nqpts, 476 interp, grad, qref, qweight) 477 478 # CeedQFunction 479 def QFunction(self, vlength, f, source): 480 """Ceed QFunction: point-wise operation at quadrature points for 481 evaluating volumetric terms. 482 483 Args: 484 vlength: vector length. Caller must ensure that number of quadrature 485 points is a multiple of vlength 486 f: ctypes function pointer to evaluate action at quadrature points 487 source: absolute path to source of QFunction, 488 "\\abs_path\\file.h:function_name 489 490 Returns: 491 qfunction: Ceed QFunction""" 492 493 return QFunction(self, vlength, f, source) 494 495 def QFunctionByName(self, name): 496 """Ceed QFunction By Name: point-wise operation at quadrature points 497 from a given gallery, for evaluating volumetric terms. 498 499 Args: 500 name: name of QFunction to use from gallery 501 502 Returns: 503 qfunction: Ceed QFunction By Name""" 504 505 return QFunctionByName(self, name) 506 507 def IdentityQFunction(self, size, inmode, outmode): 508 """Ceed Idenity QFunction: identity qfunction operation. 509 510 Args: 511 size: size of the qfunction fields 512 **inmode: CeedEvalMode for input to Ceed QFunction 513 **outmode: CeedEvalMode for output to Ceed QFunction 514 515 Returns: 516 qfunction: Ceed Identity QFunction""" 517 518 return IdentityQFunction(self, size, inmode, outmode) 519 520 def QFunctionContext(self): 521 """Ceed QFunction Context: stores Ceed QFunction user context data. 522 523 Returns: 524 userContext: Ceed QFunction Context""" 525 526 return QFunctionContext(self) 527 528 # CeedOperator 529 def Operator(self, qf, dqf=None, qdfT=None): 530 """Ceed Operator: composed FE-type operations on vectors. 531 532 Args: 533 qf: QFunction defining the action of the operator at quadrature 534 points 535 **dqf: QFunction defining the action of the Jacobian, default None 536 **dqfT: QFunction defining the action of the transpose of the 537 Jacobian, default None 538 539 Returns: 540 operator: Ceed Operator""" 541 542 return Operator(self, qf, dqf, qdfT) 543 544 def CompositeOperator(self): 545 """Composite Ceed Operator: composition of multiple CeedOperators. 546 547 Returns: 548 operator: Ceed Composite Operator""" 549 550 return CompositeOperator(self) 551 552 # Destructor 553 def __del__(self): 554 # libCEED call 555 lib.CeedDestroy(self._pointer) 556 557# ------------------------------------------------------------------------------ 558