1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights 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 #include "ceed-ref.h" 20 21 static int CeedOperatorDestroy_Ref(CeedOperator op) { 22 CeedOperator_Ref *impl = op->data; 23 int ierr; 24 25 for (CeedInt i=0; i<impl->numein+impl->numeout; i++) { 26 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChk(ierr); 27 } 28 ierr = CeedFree(&impl->evecs); CeedChk(ierr); 29 ierr = CeedFree(&impl->edata); CeedChk(ierr); 30 31 for (CeedInt i=0; i<impl->numqin+impl->numqout; i++) { 32 ierr = CeedFree(&impl->qdata_alloc[i]); CeedChk(ierr); 33 } 34 ierr = CeedFree(&impl->qdata_alloc); CeedChk(ierr); 35 ierr = CeedFree(&impl->qdata); CeedChk(ierr); 36 37 ierr = CeedFree(&impl->indata); CeedChk(ierr); 38 ierr = CeedFree(&impl->outdata); CeedChk(ierr); 39 40 ierr = CeedFree(&op->data); CeedChk(ierr); 41 return 0; 42 } 43 44 /* 45 Setup infields or outfields 46 */ 47 static int CeedOperatorSetupFields_Ref(struct CeedQFunctionField qfields[16], 48 struct CeedOperatorField ofields[16], 49 CeedVector *evecs, CeedScalar **qdata, 50 CeedScalar **qdata_alloc, CeedScalar **indata, 51 CeedInt starti, CeedInt startq, 52 CeedInt numfields, CeedInt Q) { 53 CeedInt dim, ierr, iq=startq, ncomp; 54 55 // Loop over fields 56 for (CeedInt i=0; i<numfields; i++) { 57 CeedEvalMode emode = qfields[i].emode; 58 59 if (emode != CEED_EVAL_WEIGHT) { 60 ierr = CeedElemRestrictionCreateVector(ofields[i].Erestrict, NULL, 61 &evecs[i+starti]); 62 CeedChk(ierr); 63 } 64 65 switch(emode) { 66 case CEED_EVAL_NONE: 67 break; // No action 68 case CEED_EVAL_INTERP: 69 ncomp = qfields[i].ncomp; 70 ierr = CeedMalloc(Q*ncomp, &qdata_alloc[iq]); CeedChk(ierr); 71 qdata[i + starti] = qdata_alloc[iq]; 72 iq++; 73 break; 74 case CEED_EVAL_GRAD: 75 ncomp = qfields[i].ncomp; 76 dim = ofields[i].basis->dim; 77 ierr = CeedMalloc(Q*ncomp*dim, &qdata_alloc[iq]); CeedChk(ierr); 78 qdata[i + starti] = qdata_alloc[iq]; 79 iq++; 80 break; 81 case CEED_EVAL_WEIGHT: // Only on input fields 82 ierr = CeedMalloc(Q, &qdata_alloc[iq]); CeedChk(ierr); 83 ierr = CeedBasisApply(ofields[iq].basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, 84 NULL, qdata_alloc[iq]); CeedChk(ierr); 85 qdata[i] = qdata_alloc[iq]; 86 indata[i] = qdata[i]; 87 iq++; 88 break; 89 case CEED_EVAL_DIV: 90 break; // Not implimented 91 case CEED_EVAL_CURL: 92 break; // Not implimented 93 } 94 } 95 return 0; 96 } 97 98 /* 99 CeedOperator needs to connect all the named fields (be they active or passive) 100 to the named inputs and outputs of its CeedQFunction. 101 */ 102 static int CeedOperatorSetup_Ref(CeedOperator op) { 103 if (op->setupdone) return 0; 104 CeedOperator_Ref *impl = op->data; 105 CeedQFunction qf = op->qf; 106 CeedInt Q = op->numqpoints, numinputfields, numoutputfields; 107 int ierr; 108 109 // Count infield and outfield array sizes and evectors 110 ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); CeedChk(ierr); 111 impl->numein = numinputfields; 112 for (CeedInt i=0; i<numinputfields; i++) { 113 CeedEvalMode emode = qf->inputfields[i].emode; 114 impl->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) + 115 !!(emode & CEED_EVAL_WEIGHT); 116 } 117 impl->numeout = numoutputfields; 118 for (CeedInt i=0; i<numoutputfields; i++) { 119 CeedEvalMode emode = qf->outputfields[i].emode; 120 impl->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD); 121 } 122 123 // Allocate 124 ierr = CeedCalloc(impl->numein + impl->numeout, &impl->evecs); CeedChk(ierr); 125 ierr = CeedCalloc(impl->numein + impl->numeout, &impl->edata); 126 CeedChk(ierr); 127 128 ierr = CeedCalloc(impl->numqin + impl->numqout, &impl->qdata_alloc); 129 CeedChk(ierr); 130 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->qdata); 131 CeedChk(ierr); 132 133 ierr = CeedCalloc(16, &impl->indata); CeedChk(ierr); 134 ierr = CeedCalloc(16, &impl->outdata); CeedChk(ierr); 135 136 // Set up infield and outfield pointer arrays 137 // Infields 138 ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields, 139 impl->evecs, impl->qdata, impl->qdata_alloc, 140 impl->indata, 0, 0, 141 numinputfields, Q); CeedChk(ierr); 142 143 // Outfields 144 ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields, 145 impl->evecs, impl->qdata, impl->qdata_alloc, 146 impl->indata, numinputfields, 147 impl->numqin, numoutputfields, Q); CeedChk(ierr); 148 149 // Input Qvecs 150 for (CeedInt i=0; i<numinputfields; i++) { 151 CeedEvalMode emode = qf->inputfields[i].emode; 152 if ((emode != CEED_EVAL_NONE) && (emode != CEED_EVAL_WEIGHT)) 153 impl->indata[i] = impl->qdata[i]; 154 } 155 // Output Qvecs 156 for (CeedInt i=0; i<numoutputfields; i++) { 157 CeedEvalMode emode = qf->outputfields[i].emode; 158 if (emode != CEED_EVAL_NONE) 159 impl->outdata[i] = impl->qdata[i + numinputfields]; 160 } 161 162 op->setupdone = 1; 163 164 return 0; 165 } 166 167 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec, 168 CeedVector outvec, CeedRequest *request) { 169 CeedOperator_Ref *impl = op->data; 170 CeedInt Q = op->numqpoints, elemsize, numinputfields, numoutputfields; 171 int ierr; 172 CeedQFunction qf = op->qf; 173 CeedTransposeMode lmode = CEED_NOTRANSPOSE; 174 175 // Setup 176 ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr); 177 ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); CeedChk(ierr); 178 179 // Input Evecs and Restriction 180 for (CeedInt i=0; i<numinputfields; i++) { 181 CeedEvalMode emode = qf->inputfields[i].emode; 182 if (emode == CEED_EVAL_WEIGHT) { // Skip 183 } else { 184 // Active 185 if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) { 186 // Restrict 187 ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE, 188 lmode, invec, impl->evecs[i], 189 request); CeedChk(ierr); 190 // Get evec 191 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST, 192 (const CeedScalar **) &impl->edata[i]); CeedChk(ierr); 193 } else { 194 // Passive 195 // Restrict 196 ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE, 197 lmode, op->inputfields[i].vec, impl->evecs[i], 198 request); CeedChk(ierr); 199 // Get evec 200 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST, 201 (const CeedScalar **) &impl->edata[i]); CeedChk(ierr); 202 } 203 } 204 } 205 206 // Output Evecs 207 for (CeedInt i=0; i<numoutputfields; i++) { 208 ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST, 209 &impl->edata[i + numinputfields]); CeedChk(ierr); 210 } 211 212 // Loop through elements 213 for (CeedInt e=0; e<op->numelements; e++) { 214 // Input basis apply if needed 215 for (CeedInt i=0; i<numinputfields; i++) { 216 // Get elemsize, emode, ncomp 217 elemsize = op->inputfields[i].Erestrict->elemsize; 218 CeedEvalMode emode = qf->inputfields[i].emode; 219 CeedInt ncomp = qf->inputfields[i].ncomp; 220 // Basis action 221 switch(emode) { 222 case CEED_EVAL_NONE: 223 impl->indata[i] = &impl->edata[i][e*Q*ncomp]; 224 break; 225 case CEED_EVAL_INTERP: 226 ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE, 227 CEED_EVAL_INTERP, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]); 228 CeedChk(ierr); 229 break; 230 case CEED_EVAL_GRAD: 231 ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE, 232 CEED_EVAL_GRAD, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]); 233 CeedChk(ierr); 234 break; 235 case CEED_EVAL_WEIGHT: 236 break; // No action 237 case CEED_EVAL_DIV: 238 break; // Not implimented 239 case CEED_EVAL_CURL: 240 break; // Not implimented 241 } 242 } 243 // Output pointers 244 for (CeedInt i=0; i<numoutputfields; i++) { 245 CeedEvalMode emode = qf->outputfields[i].emode; 246 if (emode == CEED_EVAL_NONE) { 247 CeedInt ncomp = qf->outputfields[i].ncomp; 248 impl->outdata[i] = &impl->edata[i + numinputfields][e*Q*ncomp]; 249 } 250 } 251 // Q function 252 ierr = CeedQFunctionApply(op->qf, Q, (const CeedScalar * const*) impl->indata, 253 impl->outdata); CeedChk(ierr); 254 255 // Output basis apply if needed 256 for (CeedInt i=0; i<numoutputfields; i++) { 257 // Get elemsize, emode, ncomp 258 elemsize = op->outputfields[i].Erestrict->elemsize; 259 CeedInt ncomp = qf->outputfields[i].ncomp; 260 CeedEvalMode emode = qf->outputfields[i].emode; 261 // Basis action 262 switch(emode) { 263 case CEED_EVAL_NONE: 264 break; // No action 265 case CEED_EVAL_INTERP: 266 ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE, 267 CEED_EVAL_INTERP, impl->outdata[i], 268 &impl->edata[i + numinputfields][e*elemsize*ncomp]); 269 CeedChk(ierr); 270 break; 271 case CEED_EVAL_GRAD: 272 ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE, 273 CEED_EVAL_GRAD, 274 impl->outdata[i], &impl->edata[i + numinputfields][e*elemsize*ncomp]); 275 CeedChk(ierr); 276 break; 277 case CEED_EVAL_WEIGHT: 278 return CeedError(op->ceed, 1, 279 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 280 break; // Should not occur 281 case CEED_EVAL_DIV: 282 break; // Not implimented 283 case CEED_EVAL_CURL: 284 break; // Not implimented 285 } 286 } 287 } 288 289 // Output restriction 290 for (CeedInt i=0; i<numoutputfields; i++) { 291 // Restore evec 292 ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein], 293 &impl->edata[i + numinputfields]); CeedChk(ierr); 294 // Active 295 if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) { 296 // Zero lvec 297 ierr = CeedVectorSetValue(outvec, 0.0); CeedChk(ierr); 298 // Restrict 299 ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE, 300 lmode, impl->evecs[i+impl->numein], outvec, request); CeedChk(ierr); 301 } else { 302 // Passive 303 // Zero lvec 304 ierr = CeedVectorSetValue(op->outputfields[i].vec, 0.0); CeedChk(ierr); 305 // Restrict 306 ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE, 307 lmode, impl->evecs[i+impl->numein], op->outputfields[i].vec, 308 request); CeedChk(ierr); 309 } 310 } 311 312 // Restore input arrays 313 for (CeedInt i=0; i<numinputfields; i++) { 314 CeedEvalMode emode = qf->inputfields[i].emode; 315 if (emode == CEED_EVAL_WEIGHT) { // Skip 316 } else { 317 ierr = CeedVectorRestoreArrayRead(impl->evecs[i], 318 (const CeedScalar **) &impl->edata[i]); CeedChk(ierr); 319 } 320 } 321 322 return 0; 323 } 324 325 int CeedOperatorCreate_Ref(CeedOperator op) { 326 CeedOperator_Ref *impl; 327 int ierr; 328 329 ierr = CeedCalloc(1, &impl); CeedChk(ierr); 330 op->data = impl; 331 op->Destroy = CeedOperatorDestroy_Ref; 332 op->Apply = CeedOperatorApply_Ref; 333 return 0; 334 } 335