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