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