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 17 #include <ceed-impl.h> 18 #include <string.h> 19 20 /// @file 21 /// Implementation of public CeedOperator interfaces 22 /// 23 /// @addtogroup CeedOperator 24 /// @{ 25 26 /** 27 @brief Create an operator from element restriction, basis, and QFunction 28 29 @param ceed A Ceed object where the CeedOperator will be created 30 @param qf QFunction defining the action of the operator at quadrature points 31 @param dqf QFunction defining the action of the Jacobian of @a qf (or NULL) 32 @param dqfT QFunction defining the action of the transpose of the Jacobian 33 of @a qf (or NULL) 34 @param[out] op Address of the variable where the newly created 35 CeedOperator will be stored 36 37 @return An error code: 0 - success, otherwise - failure 38 39 @ref Basic 40 */ 41 int CeedOperatorCreate(Ceed ceed, CeedQFunction qf, CeedQFunction dqf, 42 CeedQFunction dqfT, CeedOperator *op) { 43 int ierr; 44 45 if (!ceed->OperatorCreate) return CeedError(ceed, 1, 46 "Backend does not support OperatorCreate"); 47 ierr = CeedCalloc(1,op); CeedChk(ierr); 48 (*op)->ceed = ceed; 49 ceed->refcount++; 50 (*op)->refcount = 1; 51 (*op)->qf = qf; 52 qf->refcount++; 53 (*op)->dqf = dqf; 54 if (dqf) dqf->refcount++; 55 (*op)->dqfT = dqfT; 56 if (dqfT) dqfT->refcount++; 57 ierr = ceed->OperatorCreate(*op); CeedChk(ierr); 58 return 0; 59 } 60 61 /** 62 @brief Provide a field to a CeedOperator for use by its CeedQFunction 63 64 This function is used to specify both active and passive fields to a 65 CeedOperator. For passive fields, a vector @arg v must be provided. Passive 66 fields can inputs or outputs (updated in-place when operator is applied). 67 68 Active fields must be specified using this function, but their data (in a 69 CeedVector) is passed in CeedOperatorApply(). There can be at most one active 70 input and at most one active output. 71 72 @param op Ceedoperator on which to provide the field 73 @param fieldname Name of the field (to be matched with the name used by CeedQFunction) 74 @param r CeedElemRestriction 75 @param b CeedBasis in which the field resides or CEED_BASIS_COLLOCATED 76 if collocated with quadrature points 77 @param v CeedVector to be used by CeedOperator or CEED_VECTOR_ACTIVE 78 if field is active or CEED_VECTOR_NONE if using 79 CEED_EVAL_WEIGHT in the qfunction 80 81 @return An error code: 0 - success, otherwise - failure 82 83 @ref Basic 84 **/ 85 int CeedOperatorSetField(CeedOperator op, const char *fieldname, 86 CeedElemRestriction r, CeedBasis b, 87 CeedVector v) { 88 int ierr; 89 CeedInt numelements; 90 ierr = CeedElemRestrictionGetNumElements(r, &numelements); CeedChk(ierr); 91 if (op->numelements && op->numelements != numelements) 92 return CeedError(op->ceed, 1, 93 "ElemRestriction with %d elements incompatible with prior %d elements", 94 numelements, op->numelements); 95 op->numelements = numelements; 96 97 if (b != CEED_BASIS_COLLOCATED) { 98 CeedInt numqpoints; 99 ierr = CeedBasisGetNumQuadraturePoints(b, &numqpoints); CeedChk(ierr); 100 if (op->numqpoints && op->numqpoints != numqpoints) 101 return CeedError(op->ceed, 1, 102 "Basis with %d quadrature points incompatible with prior %d points", 103 numqpoints, op->numqpoints); 104 op->numqpoints = numqpoints; 105 } 106 struct CeedOperatorField *ofield; 107 for (CeedInt i=0; i<op->qf->numinputfields; i++) { 108 if (!strcmp(fieldname, op->qf->inputfields[i].fieldname)) { 109 ofield = &op->inputfields[i]; 110 goto found; 111 } 112 } 113 for (CeedInt i=0; i<op->qf->numoutputfields; i++) { 114 if (!strcmp(fieldname, op->qf->outputfields[i].fieldname)) { 115 ofield = &op->outputfields[i]; 116 goto found; 117 } 118 } 119 return CeedError(op->ceed, 1, "QFunction has no knowledge of field '%s'", 120 fieldname); 121 found: 122 ofield->Erestrict = r; 123 ofield->basis = b; 124 ofield->vec = v; 125 op->nfields += 1; 126 return 0; 127 } 128 129 /** 130 @brief Apply CeedOperator to a vector 131 132 This computes the action of the operator on the specified (active) input, 133 yielding its (active) output. All inputs and outputs must be specified using 134 CeedOperatorSetField(). 135 136 @param op CeedOperator to apply 137 @param[in] in CeedVector containing input state or NULL if there are no 138 active inputs 139 @param[out] out CeedVector to store result of applying operator (must be 140 distinct from @a in) or NULL if there are no active outputs 141 @param request Address of CeedRequest for non-blocking completion, else 142 CEED_REQUEST_IMMEDIATE 143 144 @return An error code: 0 - success, otherwise - failure 145 146 @ref Basic 147 **/ 148 int CeedOperatorApply(CeedOperator op, CeedVector in, 149 CeedVector out, CeedRequest *request) { 150 int ierr; 151 Ceed ceed = op->ceed; 152 CeedQFunction qf = op->qf; 153 154 if (op->nfields == 0) return CeedError(ceed, 1, "No operator fields set"); 155 if (op->nfields < qf->numinputfields + qf->numoutputfields) return CeedError( 156 ceed, 1, "Not all operator fields set"); 157 if (op->numelements == 0) return CeedError(ceed, 1, 158 "At least one restriction required"); 159 if (op->numqpoints == 0) return CeedError(ceed, 1, 160 "At least one non-collocated basis required"); 161 ierr = op->Apply(op, in, out, request); CeedChk(ierr); 162 return 0; 163 } 164 165 /** 166 @brief Destroy a CeedOperator 167 168 @param op CeedOperator to destroy 169 170 @return An error code: 0 - success, otherwise - failure 171 172 @ref Basic 173 **/ 174 int CeedOperatorDestroy(CeedOperator *op) { 175 int ierr; 176 177 if (!*op || --(*op)->refcount > 0) return 0; 178 if ((*op)->Destroy) { 179 ierr = (*op)->Destroy(*op); CeedChk(ierr); 180 } 181 ierr = CeedQFunctionDestroy(&(*op)->qf); CeedChk(ierr); 182 ierr = CeedQFunctionDestroy(&(*op)->dqf); CeedChk(ierr); 183 ierr = CeedQFunctionDestroy(&(*op)->dqfT); CeedChk(ierr); 184 ierr = CeedDestroy(&(*op)->ceed); CeedChk(ierr); 185 ierr = CeedFree(op); CeedChk(ierr); 186 return 0; 187 } 188 189 /// @} 190