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 tempfile 22from abc import ABC 23from .ceed_vector import Vector 24from .ceed_basis import BasisTensorH1, BasisTensorH1Lagrange, BasisH1 25from .ceed_elemrestriction import ElemRestriction, StridedElemRestriction, BlockedElemRestriction, BlockedStridedElemRestriction 26from .ceed_qfunction import QFunction, QFunctionByName, IdentityQFunction 27from .ceed_qfunctioncontext import QFunctionContext 28from .ceed_operator import Operator, CompositeOperator 29from .ceed_constants import * 30 31# ------------------------------------------------------------------------------ 32 33 34class Ceed(): 35 """Ceed: core components.""" 36 # Attributes 37 _pointer = ffi.NULL 38 39 # Constructor 40 def __init__(self, resource="/cpu/self", on_error="store"): 41 # libCEED object 42 self._pointer = ffi.new("Ceed *") 43 44 # libCEED call 45 resourceAscii = ffi.new("char[]", resource.encode("ascii")) 46 os.environ["CEED_ERROR_HANDLER"] = "return" 47 err_code = lib.CeedInit(resourceAscii, self._pointer) 48 if err_code: 49 raise Exception("Error initializing backend resource: " + resource) 50 error_handlers = dict( 51 store="CeedErrorStore", 52 abort="CeedErrorAbort", 53 ) 54 lib.CeedSetErrorHandler( 55 self._pointer[0], ffi.addressof( 56 lib, error_handlers[on_error])) 57 58 # Representation 59 def __repr__(self): 60 return "<Ceed instance at " + hex(id(self)) + ">" 61 62 # String conversion for print() to stdout 63 def __str__(self): 64 """View a Ceed via print().""" 65 66 # libCEED call 67 with tempfile.NamedTemporaryFile() as key_file: 68 with open(key_file.name, 'r+') as stream_file: 69 stream = ffi.cast("FILE *", stream_file) 70 71 err_code = lib.CeedView(self._pointer[0], stream) 72 self._check_error(err_code) 73 74 stream_file.seek(0) 75 out_string = stream_file.read() 76 77 return out_string 78 79 # Error handler 80 def _check_error(self, err_code): 81 """Check return code and retrieve error message for non-zero code""" 82 if (err_code): 83 message = ffi.new("char **") 84 lib.CeedGetErrorMessage(self._pointer[0], message) 85 raise Exception(ffi.string(message[0]).decode("UTF-8")) 86 87 # Get Resource 88 def get_resource(self): 89 """Get the full resource name for a Ceed context. 90 91 Returns: 92 resource: resource name""" 93 94 # libCEED call 95 resource = ffi.new("char **") 96 err_code = lib.CeedGetResource(self._pointer[0], resource) 97 self._check_error(err_code) 98 99 return ffi.string(resource[0]).decode("UTF-8") 100 101 # Preferred MemType 102 def get_preferred_memtype(self): 103 """Return Ceed preferred memory type. 104 105 Returns: 106 memtype: Ceed preferred memory type""" 107 108 # libCEED call 109 memtype = ffi.new("CeedMemType *", MEM_HOST) 110 err_code = lib.CeedGetPreferredMemType(self._pointer[0], memtype) 111 self._check_error(err_code) 112 113 return memtype[0] 114 115 # CeedVector 116 def Vector(self, size): 117 """Ceed Vector: storing and manipulating vectors. 118 119 Args: 120 size: length of vector 121 122 Returns: 123 vector: Ceed Vector""" 124 125 return Vector(self, size) 126 127 # CeedElemRestriction 128 def ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, offsets, 129 memtype=lib.CEED_MEM_HOST, cmode=lib.CEED_COPY_VALUES): 130 """Ceed ElemRestriction: restriction from local vectors to elements. 131 132 Args: 133 nelem: number of elements described by the restriction 134 elemsize: size (number of nodes) per element 135 ncomp: number of field components per interpolation node 136 (1 for scalar fields) 137 compstride: Stride between components for the same L-vector "node". 138 Data for node i, component k can be found in the 139 L-vector at index [offsets[i] + k*compstride]. 140 lsize: The size of the L-vector. This vector may be larger than 141 the elements and fields given by this restriction. 142 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 143 ordered list of the offsets (into the input Ceed Vector) 144 for the unknowns corresponding to element i, where 145 0 <= i < nelem. All offsets must be in the range 146 [0, lsize - 1]. 147 **memtype: memory type of the offsets array, default CEED_MEM_HOST 148 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 149 150 Returns: 151 elemrestriction: Ceed ElemRestiction""" 152 153 return ElemRestriction(self, nelem, elemsize, ncomp, compstride, lsize, 154 offsets, memtype=memtype, cmode=cmode) 155 156 def StridedElemRestriction(self, nelem, elemsize, ncomp, lsize, strides): 157 """Ceed Identity ElemRestriction: strided restriction from local vectors 158 to elements. 159 160 Args: 161 nelem: number of elements described by the restriction 162 elemsize: size (number of nodes) per element 163 ncomp: number of field components per interpolation node 164 (1 for scalar fields) 165 lsize: The size of the L-vector. This vector may be larger than 166 the elements and fields given by this restriction. 167 *strides: Array for strides between [nodes, components, elements]. 168 The data for node i, component j, element k in the 169 L-vector is given by 170 i*strides[0] + j*strides[1] + k*strides[2] 171 172 Returns: 173 elemrestriction: Ceed Strided ElemRestiction""" 174 175 return StridedElemRestriction( 176 self, nelem, elemsize, ncomp, lsize, strides) 177 178 def BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, compstride, 179 lsize, offsets, memtype=lib.CEED_MEM_HOST, 180 cmode=lib.CEED_COPY_VALUES): 181 """Ceed Blocked ElemRestriction: blocked restriction from local vectors 182 to elements. 183 184 Args: 185 nelem: number of elements described by the restriction 186 elemsize: size (number of nodes) per element 187 blksize: number of elements in a block 188 ncomp: number of field components per interpolation node 189 (1 for scalar fields) 190 lsize: The size of the L-vector. This vector may be larger than 191 the elements and fields given by this restriction. 192 *offsets: Numpy array of shape [nelem, elemsize]. Row i holds the 193 ordered list of the offsets (into the input Ceed Vector) 194 for the unknowns corresponding to element i, where 195 0 <= i < nelem. All offsets must be in the range 196 [0, lsize - 1]. The backend will permute and pad this 197 array to the desired ordering for the blocksize, which is 198 typically given by the backend. The default reordering is 199 to interlace elements. 200 **memtype: memory type of the offsets array, default CEED_MEM_HOST 201 **cmode: copy mode for the offsets array, default CEED_COPY_VALUES 202 203 Returns: 204 elemrestriction: Ceed Blocked ElemRestiction""" 205 206 return BlockedElemRestriction(self, nelem, elemsize, blksize, ncomp, 207 compstride, lsize, offsets, 208 memtype=memtype, cmode=cmode) 209 210 def BlockedStridedElemRestriction(self, nelem, elemsize, blksize, ncomp, 211 lsize, strides): 212 """Ceed Blocked Strided ElemRestriction: blocked and strided restriction 213 from local vectors to elements. 214 215 Args: 216 nelem: number of elements described in the restriction 217 elemsize: size (number of nodes) per element 218 blksize: number of elements in a block 219 ncomp: number of field components per interpolation node 220 (1 for scalar fields) 221 lsize: The size of the L-vector. This vector may be larger than 222 the elements and fields given by this restriction. 223 *strides: Array for strides between [nodes, components, elements]. 224 The data for node i, component j, element k in the 225 L-vector is given by 226 i*strides[0] + j*strides[1] + k*strides[2] 227 228 Returns: 229 elemrestriction: Ceed Strided ElemRestiction""" 230 231 return BlockedStridedElemRestriction(self, nelem, elemsize, blksize, 232 ncomp, lsize, strides) 233 234 # CeedBasis 235 def BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 236 qref1d, qweight1d): 237 """Ceed Tensor H1 Basis: finite element tensor-product basis objects for 238 H^1 discretizations. 239 240 Args: 241 dim: topological dimension 242 ncomp: number of field components (1 for scalar fields) 243 P1d: number of nodes in one dimension 244 Q1d: number of quadrature points in one dimension 245 *interp1d: Numpy array holding the row-major (Q1d * P1d) matrix 246 expressing the values of nodal basis functions at 247 quadrature points 248 *grad1d: Numpy array holding the row-major (Q1d * P1d) matrix 249 expressing the derivatives of nodal basis functions at 250 quadrature points 251 *qref1d: Numpy array of length Q1d holding the locations of 252 quadrature points on the 1D reference element [-1, 1] 253 *qweight1d: Numpy array of length Q1d holding the quadrature 254 weights on the reference element 255 256 Returns: 257 basis: Ceed Basis""" 258 259 return BasisTensorH1(self, dim, ncomp, P1d, Q1d, interp1d, grad1d, 260 qref1d, qweight1d) 261 262 def BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode): 263 """Ceed Tensor H1 Lagrange Basis: finite element tensor-product Lagrange 264 basis objects for H^1 discretizations. 265 266 Args: 267 dim: topological dimension 268 ncomp: number of field components (1 for scalar fields) 269 P: number of Gauss-Lobatto nodes in one dimension. The 270 polynomial degree of the resulting Q_k element is k=P-1. 271 Q: number of quadrature points in one dimension 272 qmode: distribution of the Q quadrature points (affects order of 273 accuracy for the quadrature) 274 275 Returns: 276 basis: Ceed Basis""" 277 278 return BasisTensorH1Lagrange(self, dim, ncomp, P, Q, qmode) 279 280 def BasisH1(self, topo, ncomp, nnodes, nqpts, interp, grad, qref, qweight): 281 """Ceed H1 Basis: finite element non tensor-product basis for H^1 282 discretizations. 283 284 Args: 285 topo: topology of the element, e.g. hypercube, simplex, etc 286 ncomp: number of field components (1 for scalar fields) 287 nnodes: total number of nodes 288 nqpts: total number of quadrature points 289 *interp: Numpy array holding the row-major (nqpts * nnodes) matrix 290 expressing the values of nodal basis functions at 291 quadrature points 292 *grad: Numpy array holding the row-major (nqpts * dim * nnodes) 293 matrix expressing the derivatives of nodal basis functions 294 at quadrature points 295 *qref: Numpy array of length (nqpts * dim) holding the locations of 296 quadrature points on the reference element [-1, 1] 297 *qweight: Numpy array of length nnodes holding the quadrature 298 weights on the reference element 299 300 Returns: 301 basis: Ceed Basis""" 302 303 return BasisH1(self, topo, ncomp, nnodes, nqpts, 304 interp, grad, qref, qweight) 305 306 # CeedQFunction 307 def QFunction(self, vlength, f, source): 308 """Ceed QFunction: point-wise operation at quadrature points for 309 evaluating volumetric terms. 310 311 Args: 312 vlength: vector length. Caller must ensure that number of quadrature 313 points is a multiple of vlength 314 f: ctypes function pointer to evaluate action at quadrature points 315 source: absolute path to source of QFunction, 316 "\\abs_path\\file.h:function_name 317 318 Returns: 319 qfunction: Ceed QFunction""" 320 321 return QFunction(self, vlength, f, source) 322 323 def QFunctionByName(self, name): 324 """Ceed QFunction By Name: point-wise operation at quadrature points 325 from a given gallery, for evaluating volumetric terms. 326 327 Args: 328 name: name of QFunction to use from gallery 329 330 Returns: 331 qfunction: Ceed QFunction By Name""" 332 333 return QFunctionByName(self, name) 334 335 def IdentityQFunction(self, size, inmode, outmode): 336 """Ceed Idenity QFunction: identity qfunction operation. 337 338 Args: 339 size: size of the qfunction fields 340 **inmode: CeedEvalMode for input to Ceed QFunction 341 **outmode: CeedEvalMode for output to Ceed QFunction 342 343 Returns: 344 qfunction: Ceed Identity QFunction""" 345 346 return IdentityQFunction(self, size, inmode, outmode) 347 348 def QFunctionContext(self): 349 """Ceed QFunction Context: stores Ceed QFunction user context data. 350 351 Returns: 352 userContext: Ceed QFunction Context""" 353 354 return QFunctionContext(self) 355 356 # CeedOperator 357 def Operator(self, qf, dqf=None, qdfT=None): 358 """Ceed Operator: composed FE-type operations on vectors. 359 360 Args: 361 qf: QFunction defining the action of the operator at quadrature 362 points 363 **dqf: QFunction defining the action of the Jacobian, default None 364 **dqfT: QFunction defining the action of the transpose of the 365 Jacobian, default None 366 367 Returns: 368 operator: Ceed Operator""" 369 370 return Operator(self, qf, dqf, qdfT) 371 372 def CompositeOperator(self): 373 """Composite Ceed Operator: composition of multiple CeedOperators. 374 375 Returns: 376 operator: Ceed Composite Operator""" 377 378 return CompositeOperator(self) 379 380 # Destructor 381 def __del__(self): 382 # libCEED call 383 lib.CeedDestroy(self._pointer) 384 385# ------------------------------------------------------------------------------ 386