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 tempfile 19from abc import ABC 20from .ceed_constants import REQUEST_IMMEDIATE, REQUEST_ORDERED, NOTRANSPOSE 21 22# ------------------------------------------------------------------------------ 23 24 25class _OperatorBase(ABC): 26 """Ceed Operator: composed FE-type operations on vectors.""" 27 28 # Attributes 29 _ceed = None 30 _pointer = ffi.NULL 31 32 # Destructor 33 def __del__(self): 34 # libCEED call 35 err_code = lib.CeedOperatorDestroy(self._pointer) 36 self._ceed._check_error(err_code) 37 38 # Representation 39 def __repr__(self): 40 return "<CeedOperator instance at " + hex(id(self)) + ">" 41 42 # String conversion for print() to stdout 43 def __str__(self): 44 """View an Operator via print().""" 45 46 # libCEED call 47 with tempfile.NamedTemporaryFile() as key_file: 48 with open(key_file.name, 'r+') as stream_file: 49 stream = ffi.cast("FILE *", stream_file) 50 51 err_code = lib.CeedOperatorView(self._pointer[0], stream) 52 self._ceed._check_error(err_code) 53 54 stream_file.seek(0) 55 out_string = stream_file.read() 56 57 return out_string 58 59 # Assemble linear diagonal 60 def linear_assemble_diagonal(self, d, request=REQUEST_IMMEDIATE): 61 """Assemble the diagonal of a square linear Operator 62 63 Args: 64 d: Vector to store assembled Operator diagonal 65 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 66 67 # libCEED call 68 err_code = lib.CeedOperatorLinearAssembleDiagonal(self._pointer[0], 69 d._pointer[0], request) 70 self._ceed._check_error(err_code) 71 72 # Assemble add linear diagonal 73 def linear_assemble_add_diagonal(self, d, request=REQUEST_IMMEDIATE): 74 """Sum the diagonal of a square linear Operator into a Vector 75 76 Args: 77 d: Vector to store assembled Operator diagonal 78 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 79 80 # libCEED call 81 err_code = lib.CeedOperatorLinearAssembleAddDiagonal(self._pointer[0], 82 d._pointer[0], request) 83 self._ceed._check_error(err_code) 84 85 # Assemble linear point block diagonal 86 def linear_assemble_point_block_diagonal( 87 self, d, request=REQUEST_IMMEDIATE): 88 """Assemble the point block diagonal of a square linear Operator 89 90 Args: 91 d: Vector to store assembled Operator point block diagonal, 92 provided in row-major form with an ncomp*ncomp block 93 at each node 94 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 95 96 # libCEED call 97 err_code = lib.CeedOperatorLinearAssemblePointBlockDiagonal(self._pointer[0], 98 d._pointer[0], request) 99 self._ceed._check_error(err_code) 100 101 # Assemble linear point block diagonal 102 def linear_assemble_add_point_block_diagonal( 103 self, d, request=REQUEST_IMMEDIATE): 104 """Sum the point block diagonal of a square linear Operator into a Vector 105 106 Args: 107 d: Vector to store assembled Operator point block diagonal, 108 provided in row-major form with an ncomp*ncomp block 109 at each node 110 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 111 112 # libCEED call 113 err_code = lib.CeedOperatorLinearAssembleAddPointBlockDiagonal(self._pointer[0], 114 d._pointer[0], request) 115 self._ceed._check_error(err_code) 116 117 # Apply CeedOperator 118 def apply(self, u, v, request=REQUEST_IMMEDIATE): 119 """Apply Operator to a vector. 120 121 Args: 122 u: Vector containing input state or CEED_VECTOR_NONE if there are no 123 active inputs 124 v: Vector to store result of applying operator (must be distinct from u) 125 or CEED_VECTOR_NONE if there are no active outputs 126 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 127 128 # libCEED call 129 err_code = lib.CeedOperatorApply(self._pointer[0], u._pointer[0], v._pointer[0], 130 request) 131 self._ceed._check_error(err_code) 132 133 # Apply CeedOperator 134 def apply_add(self, u, v, request=REQUEST_IMMEDIATE): 135 """Apply Operator to a vector and add result to output vector. 136 137 Args: 138 u: Vector containing input state or CEED_VECTOR_NONE if there are no 139 active inputs 140 v: Vector to sum in result of applying operator (must be distinct from u) 141 or CEED_VECTOR_NONE if there are no active outputs 142 **request: Ceed request, default CEED_REQUEST_IMMEDIATE""" 143 144 # libCEED call 145 err_code = lib.CeedOperatorApplyAdd(self._pointer[0], u._pointer[0], v._pointer[0], 146 request) 147 self._ceed._check_error(err_code) 148 149 # Create Multigrid Level 150 def multigrid_create(self, p_mult_fine, rstr_coarse, basis_coarse): 151 """ Create a multigrid coarse operator and level transfer operators 152 for a CeedOperator with a Lagrange tensor basis for the active basis 153 154 Args: 155 p_mult_fine: L-vector multiplicity in parallel gather/scatter 156 basis_coarse: Coarse grid active vector basis 157 degree_coarse: Coarse grid basis polynomial order""" 158 159 # Operator pointers 160 opCoarsePointer = ffi.new("CeedOperator *") 161 opProlongPointer = ffi.new("CeedOperator *") 162 opRestrictPointer = ffi.new("CeedOperator *") 163 164 # libCEED call 165 lib.CeedOperatorMultigridLevelCreate(self._pointer[0], 166 p_mult_fine._pointer[0], 167 rstr_coarse._pointer[0], 168 basis_coarse._pointer[0], 169 opCoarsePointer, 170 opProlongPointer, 171 opRestrictPointer) 172 173 # Wrap operators 174 opCoarse = _OperatorWrap( 175 self._ceed, opCoarsePointer) 176 opProlong = _OperatorWrap( 177 self._ceed, opProlongPointer) 178 opRestrict = _OperatorWrap( 179 self._ceed, opRestrictPointer) 180 181 # Return 182 return [opCoarse, opProlong, opRestrict] 183 184 # Create Multigrid Level 185 def multigrid_create_tensor_h1(self, p_mult_fine, rstr_coarse, basis_coarse, 186 interp_C_to_F): 187 """ Create a multigrid coarse operator and level transfer operators 188 for a CeedOperator with a non-tensor basis for the active basis 189 190 Args: 191 p_mult_fine: L-vector multiplicity in parallel gather/scatter 192 rstr_coarse: Coarse grid restriction 193 basis_coarse: Coarse grid active vector basis 194 interp_C_to_F: Matrix for coarse to fine interpolation""" 195 196 # Setup arguments 197 interpCtoF_pointer = ffi.new("CeedScalar *") 198 interpCtoF_pointer = ffi.cast( 199 "CeedScalar *", 200 interp_C_to_F.__array_interface__['data'][0]) 201 202 # Operator pointers 203 opCoarsePointer = ffi.new("CeedOperator *") 204 opProlongPointer = ffi.new("CeedOperator *") 205 opRestrictPointer = ffi.new("CeedOperator *") 206 207 # libCEED call 208 lib.CeedOperatorMultigridLevelCreateTensorH1(self._pointer[0], 209 p_mult_fine._pointer[0], 210 rstr_coarse._pointer[0], 211 basis_coarse._pointer[0], 212 interpCtoF_pointer, 213 opCoarsePointer, 214 opProlongPointer, 215 opRestrictPointer) 216 217 # Wrap operators 218 opCoarse = _OperatorWrap( 219 self._ceed, opCoarsePointer) 220 opProlong = _OperatorWrap( 221 self._ceed, opProlongPointer) 222 opRestrict = _OperatorWrap( 223 self._ceed, opRestrictPointer) 224 225 # Return 226 return [opCoarse, opProlong, opRestrict] 227 228 # Create Multigrid Level 229 def multigrid_create_h1(self, p_mult_fine, rstr_coarse, basis_coarse, 230 interp_C_to_F): 231 """ Create a multigrid coarse operator and level transfer operators 232 for a CeedOperator with a Lagrange tensor basis for the active basis 233 234 Args: 235 p_mult_fine: L-vector multiplicity in parallel gather/scatter 236 rstr_coarse: Coarse grid restriction 237 basis_coarse: Coarse grid active vector basis 238 interp_C_to_F: Matrix for coarse to fine interpolation""" 239 240 # Setup arguments 241 interpCtoF_pointer = ffi.new("CeedScalar *") 242 interpCtoF_pointer = ffi.cast( 243 "CeedScalar *", 244 interp_C_to_F.__array_interface__['data'][0]) 245 246 # Operator pointers 247 opCoarsePointer = ffi.new("CeedOperator *") 248 opProlongPointer = ffi.new("CeedOperator *") 249 opRestrictPointer = ffi.new("CeedOperator *") 250 251 # libCEED call 252 lib.CeedOperatorMultigridLevelCreateH1(self._pointer[0], 253 p_mult_fine._pointer[0], 254 rstr_coarse._pointer[0], 255 basis_coarse._pointer[0], 256 interpCtoF_pointer, 257 opCoarsePointer, 258 opProlongPointer, 259 opRestrictPointer) 260 261 # Wrap operators 262 opCoarse = _OperatorWrap( 263 self._ceed, opCoarsePointer) 264 opProlong = _OperatorWrap( 265 self._ceed, opProlongPointer) 266 opRestrict = _OperatorWrap( 267 self._ceed, opRestrictPointer) 268 269 # Return 270 return [opCoarse, opProlong, opRestrict] 271 272 273# ------------------------------------------------------------------------------ 274 275 276class Operator(_OperatorBase): 277 """Ceed Operator: composed FE-type operations on vectors.""" 278 279 # Constructor 280 def __init__(self, ceed, qf, dqf=None, dqfT=None): 281 # CeedOperator object 282 self._pointer = ffi.new("CeedOperator *") 283 284 # Reference to Ceed 285 self._ceed = ceed 286 287 # libCEED call 288 err_code = lib.CeedOperatorCreate(self._ceed._pointer[0], qf._pointer[0], 289 dqf._pointer[0] if dqf else ffi.NULL, 290 dqfT._pointer[0] if dqfT else ffi.NULL, 291 self._pointer) 292 self._ceed._check_error(err_code) 293 294 # Add field to CeedOperator 295 def set_field(self, fieldname, restriction, basis, vector): 296 """Provide a field to an Operator for use by its QFunction. 297 298 Args: 299 fieldname: name of the field (to be matched with the same name used 300 by QFunction) 301 restriction: ElemRestriction 302 basis: Basis in which the field resides or CEED_BASIS_COLLOCATED 303 if collocated with quadrature points 304 vector: Vector to be used by Operator or CEED_VECTOR_ACTIVE 305 if field is active or CEED_VECTOR_NONE if using 306 CEED_EVAL_WEIGHT in the QFunction""" 307 308 # libCEED call 309 fieldnameAscii = ffi.new("char[]", fieldname.encode('ascii')) 310 err_code = lib.CeedOperatorSetField(self._pointer[0], fieldnameAscii, 311 restriction._pointer[0], basis._pointer[0], 312 vector._pointer[0]) 313 self._ceed._check_error(err_code) 314 315# ------------------------------------------------------------------------------ 316 317 318class CompositeOperator(_OperatorBase): 319 """Ceed Composite Operator: composition of multiple Operators.""" 320 321 # Constructor 322 def __init__(self, ceed): 323 # CeedOperator object 324 self._pointer = ffi.new("CeedOperator *") 325 326 # Reference to Ceed 327 self._ceed = ceed 328 # libCEED call 329 err_code = lib.CeedCompositeOperatorCreate( 330 self._ceed._pointer[0], self._pointer) 331 self._ceed._check_error(err_code) 332 333 # Add sub operators 334 def add_sub(self, subop): 335 """Add a sub-operator to a composite CeedOperator. 336 337 Args: 338 subop: sub-operator Operator""" 339 340 # libCEED call 341 err_code = lib.CeedCompositeOperatorAddSub( 342 self._pointer[0], subop._pointer[0]) 343 self._ceed._check_error(err_code) 344 345# ------------------------------------------------------------------------------ 346 347 348class _OperatorWrap(Operator): 349 """Wrap a CeedOperator pointer in a Operator object.""" 350 351 # Constructor 352 def __init__(self, ceed, pointer): 353 # CeedOperator object 354 self._pointer = pointer 355 356 # Reference to Ceed 357 self._ceed = ceed 358 359# ------------------------------------------------------------------------------ 360