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