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