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 <string.h> 18 #include "ceed-blocked.h" 19 #include "../ref/ceed-ref.h" 20 21 static int CeedOperatorDestroy_Blocked(CeedOperator op) { 22 int ierr; 23 CeedOperator_Blocked *impl; 24 ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr); 25 26 for (CeedInt i=0; i<impl->numein+impl->numeout; i++) { 27 ierr = CeedElemRestrictionDestroy(&impl->blkrestr[i]); CeedChk(ierr); 28 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChk(ierr); 29 } 30 ierr = CeedFree(&impl->blkrestr); CeedChk(ierr); 31 ierr = CeedFree(&impl->evecs); CeedChk(ierr); 32 ierr = CeedFree(&impl->edata); CeedChk(ierr); 33 ierr = CeedFree(&impl->inputstate); CeedChk(ierr); 34 35 for (CeedInt i=0; i<impl->numein; i++) { 36 ierr = CeedVectorDestroy(&impl->evecsin[i]); CeedChk(ierr); 37 ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChk(ierr); 38 } 39 ierr = CeedFree(&impl->evecsin); CeedChk(ierr); 40 ierr = CeedFree(&impl->qvecsin); CeedChk(ierr); 41 42 for (CeedInt i=0; i<impl->numeout; i++) { 43 ierr = CeedVectorDestroy(&impl->evecsout[i]); CeedChk(ierr); 44 ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChk(ierr); 45 } 46 ierr = CeedFree(&impl->evecsout); CeedChk(ierr); 47 ierr = CeedFree(&impl->qvecsout); CeedChk(ierr); 48 49 ierr = CeedFree(&impl); CeedChk(ierr); 50 return 0; 51 } 52 53 /* 54 Setup infields or outfields 55 */ 56 static int CeedOperatorSetupFields_Blocked(CeedQFunction qf, CeedOperator op, 57 bool inOrOut, 58 CeedElemRestriction *blkrestr, 59 CeedVector *fullevecs, CeedVector *evecs, 60 CeedVector *qvecs, CeedInt starte, 61 CeedInt numfields, CeedInt Q) { 62 CeedInt dim, ierr, ncomp, P; 63 Ceed ceed; 64 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 65 CeedBasis basis; 66 CeedElemRestriction r; 67 CeedOperatorField *opfields; 68 CeedQFunctionField *qffields; 69 if (inOrOut) { 70 ierr = CeedOperatorGetFields(op, NULL, &opfields); 71 CeedChk(ierr); 72 ierr = CeedQFunctionGetFields(qf, NULL, &qffields); 73 CeedChk(ierr); 74 } else { 75 ierr = CeedOperatorGetFields(op, &opfields, NULL); 76 CeedChk(ierr); 77 ierr = CeedQFunctionGetFields(qf, &qffields, NULL); 78 CeedChk(ierr); 79 } 80 const CeedInt blksize = 8; 81 82 // Loop over fields 83 for (CeedInt i=0; i<numfields; i++) { 84 CeedEvalMode emode; 85 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChk(ierr); 86 87 if (emode != CEED_EVAL_WEIGHT) { 88 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &r); 89 CeedChk(ierr); 90 CeedElemRestriction_Ref *data; 91 ierr = CeedElemRestrictionGetData(r, (void *)&data); 92 Ceed ceed; 93 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChk(ierr); 94 CeedInt nelem, elemsize, ndof; 95 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 96 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 97 ierr = CeedElemRestrictionGetNumDoF(r, &ndof); CeedChk(ierr); 98 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); CeedChk(ierr); 99 ierr = CeedElemRestrictionCreateBlocked(ceed, nelem, elemsize, 100 blksize, ndof, ncomp, 101 CEED_MEM_HOST, CEED_COPY_VALUES, 102 data->indices, &blkrestr[i+starte]); 103 CeedChk(ierr); 104 ierr = CeedElemRestrictionCreateVector(blkrestr[i+starte], NULL, 105 &fullevecs[i+starte]); 106 CeedChk(ierr); 107 } 108 109 switch(emode) { 110 case CEED_EVAL_NONE: 111 ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp); 112 CeedChk(ierr); 113 ierr = CeedVectorCreate(ceed, Q*ncomp*blksize, &qvecs[i]); CeedChk(ierr); 114 break; 115 case CEED_EVAL_INTERP: 116 ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp); 117 CeedChk(ierr); 118 ierr = CeedElemRestrictionGetElementSize(r, &P); 119 CeedChk(ierr); 120 ierr = CeedVectorCreate(ceed, P*ncomp*blksize, &evecs[i]); CeedChk(ierr); 121 ierr = CeedVectorCreate(ceed, Q*ncomp*blksize, &qvecs[i]); CeedChk(ierr); 122 break; 123 case CEED_EVAL_GRAD: 124 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr); 125 ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp); 126 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 127 ierr = CeedElemRestrictionGetElementSize(r, &P); 128 CeedChk(ierr); 129 ierr = CeedVectorCreate(ceed, P*ncomp*blksize, &evecs[i]); CeedChk(ierr); 130 ierr = CeedVectorCreate(ceed, Q*ncomp*dim*blksize, &qvecs[i]); CeedChk(ierr); 131 break; 132 case CEED_EVAL_WEIGHT: // Only on input fields 133 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr); 134 ierr = CeedVectorCreate(ceed, Q*blksize, &qvecs[i]); CeedChk(ierr); 135 ierr = CeedBasisApply(basis, blksize, CEED_NOTRANSPOSE, 136 CEED_EVAL_WEIGHT, NULL, qvecs[i]); CeedChk(ierr); 137 138 break; 139 case CEED_EVAL_DIV: 140 break; // Not implimented 141 case CEED_EVAL_CURL: 142 break; // Not implimented 143 } 144 } 145 return 0; 146 } 147 148 /* 149 CeedOperator needs to connect all the named fields (be they active or passive) 150 to the named inputs and outputs of its CeedQFunction. 151 */ 152 static int CeedOperatorSetup_Blocked(CeedOperator op) { 153 int ierr; 154 bool setupdone; 155 ierr = CeedOperatorGetSetupStatus(op, &setupdone); CeedChk(ierr); 156 if (setupdone) return 0; 157 Ceed ceed; 158 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 159 CeedOperator_Blocked *impl; 160 ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr); 161 CeedQFunction qf; 162 ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); 163 CeedInt Q, numinputfields, numoutputfields; 164 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr); 165 ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 166 CeedChk(ierr); 167 CeedOperatorField *opinputfields, *opoutputfields; 168 ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields); 169 CeedChk(ierr); 170 CeedQFunctionField *qfinputfields, *qfoutputfields; 171 ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields); 172 CeedChk(ierr); 173 174 // Allocate 175 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->blkrestr); 176 CeedChk(ierr); 177 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs); 178 CeedChk(ierr); 179 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->edata); 180 CeedChk(ierr); 181 182 ierr = CeedCalloc(16, &impl->inputstate); CeedChk(ierr); 183 ierr = CeedCalloc(16, &impl->evecsin); CeedChk(ierr); 184 ierr = CeedCalloc(16, &impl->evecsout); CeedChk(ierr); 185 ierr = CeedCalloc(16, &impl->qvecsin); CeedChk(ierr); 186 ierr = CeedCalloc(16, &impl->qvecsout); CeedChk(ierr); 187 188 impl->numein = numinputfields; impl->numeout = numoutputfields; 189 190 // Set up infield and outfield pointer arrays 191 // Infields 192 ierr = CeedOperatorSetupFields_Blocked(qf, op, 0, impl->blkrestr, 193 impl->evecs, impl->evecsin, 194 impl->qvecsin, 0, 195 numinputfields, Q); 196 CeedChk(ierr); 197 // Outfields 198 ierr = CeedOperatorSetupFields_Blocked(qf, op, 1, impl->blkrestr, 199 impl->evecs, impl->evecsout, 200 impl->qvecsout, numinputfields, 201 numoutputfields, Q); 202 CeedChk(ierr); 203 204 ierr = CeedOperatorSetSetupDone(op); CeedChk(ierr); 205 206 return 0; 207 } 208 209 static int CeedOperatorApply_Blocked(CeedOperator op, CeedVector invec, 210 CeedVector outvec, CeedRequest *request) { 211 int ierr; 212 CeedOperator_Blocked *impl; 213 ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr); 214 const CeedInt blksize = 8; 215 CeedInt Q, elemsize, numinputfields, numoutputfields, numelements, ncomp; 216 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr); 217 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr); 218 CeedInt nblks = (numelements/blksize) + !!(numelements%blksize); 219 CeedQFunction qf; 220 ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr); 221 ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 222 CeedChk(ierr); 223 CeedTransposeMode lmode; 224 CeedOperatorField *opinputfields, *opoutputfields; 225 ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields); 226 CeedChk(ierr); 227 CeedQFunctionField *qfinputfields, *qfoutputfields; 228 ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields); 229 CeedChk(ierr); 230 CeedEvalMode emode; 231 CeedVector vec; 232 CeedBasis basis; 233 CeedElemRestriction Erestrict; 234 uint64_t state; 235 236 // Setup 237 ierr = CeedOperatorSetup_Blocked(op); CeedChk(ierr); 238 239 // Input Evecs and Restriction 240 for (CeedInt i=0; i<numinputfields; i++) { 241 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 242 CeedChk(ierr); 243 if (emode == CEED_EVAL_WEIGHT) { // Skip 244 } else { 245 // Get input vector 246 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr); 247 if (vec == CEED_VECTOR_ACTIVE) 248 vec = invec; 249 // Restrict 250 ierr = CeedVectorGetState(vec, &state); CeedChk(ierr); 251 if (state != impl->inputstate[i] || vec == invec) { 252 ierr = CeedOperatorFieldGetLMode(opinputfields[i], &lmode); CeedChk(ierr); 253 ierr = CeedElemRestrictionApply(impl->blkrestr[i], CEED_NOTRANSPOSE, 254 lmode, vec, impl->evecs[i], 255 request); CeedChk(ierr); CeedChk(ierr); 256 impl->inputstate[i] = state; 257 } 258 // Get evec 259 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST, 260 (const CeedScalar **) &impl->edata[i]); 261 CeedChk(ierr); 262 } 263 } 264 265 // Output Evecs 266 for (CeedInt i=0; i<numoutputfields; i++) { 267 ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST, 268 &impl->edata[i + numinputfields]); CeedChk(ierr); 269 } 270 271 // Loop through elements 272 for (CeedInt e=0; e<nblks*blksize; e+=blksize) { 273 // Input basis apply if needed 274 for (CeedInt i=0; i<numinputfields; i++) { 275 // Get elemsize, emode, ncomp 276 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 277 CeedChk(ierr); 278 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 279 CeedChk(ierr); 280 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 281 CeedChk(ierr); 282 ierr = CeedQFunctionFieldGetNumComponents(qfinputfields[i], &ncomp); 283 CeedChk(ierr); 284 // Basis action 285 switch(emode) { 286 case CEED_EVAL_NONE: 287 ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_HOST, 288 CEED_USE_POINTER, 289 &impl->edata[i][e*Q*ncomp]); CeedChk(ierr); 290 break; 291 case CEED_EVAL_INTERP: 292 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr); 293 ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST, 294 CEED_USE_POINTER, 295 &impl->edata[i][e*elemsize*ncomp]); 296 CeedChk(ierr); 297 ierr = CeedBasisApply(basis, blksize, CEED_NOTRANSPOSE, 298 CEED_EVAL_INTERP, impl->evecsin[i], 299 impl->qvecsin[i]); CeedChk(ierr); 300 break; 301 case CEED_EVAL_GRAD: 302 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr); 303 ierr = CeedVectorSetArray(impl->evecsin[i], CEED_MEM_HOST, 304 CEED_USE_POINTER, 305 &impl->edata[i][e*elemsize*ncomp]); 306 CeedChk(ierr); 307 ierr = CeedBasisApply(basis, blksize, CEED_NOTRANSPOSE, 308 CEED_EVAL_GRAD, impl->evecsin[i], 309 impl->qvecsin[i]); CeedChk(ierr); 310 break; 311 case CEED_EVAL_WEIGHT: 312 break; // No action 313 case CEED_EVAL_DIV: 314 break; // Not implimented 315 case CEED_EVAL_CURL: 316 break; // Not implimented 317 } 318 } 319 320 // Output pointers 321 for (CeedInt i=0; i<numoutputfields; i++) { 322 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 323 CeedChk(ierr); 324 if (emode == CEED_EVAL_NONE) { 325 ierr = CeedQFunctionFieldGetNumComponents(qfoutputfields[i], &ncomp); 326 CeedChk(ierr); 327 ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_HOST, 328 CEED_USE_POINTER, 329 &impl->edata[i + numinputfields][e*Q*ncomp]); 330 CeedChk(ierr); 331 } 332 } 333 // Q function 334 ierr = CeedQFunctionApply(qf, Q*blksize, impl->qvecsin, impl->qvecsout); 335 CeedChk(ierr); 336 337 // Output basis apply if needed 338 for (CeedInt i=0; i<numoutputfields; i++) { 339 // Get elemsize, emode, ncomp 340 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 341 CeedChk(ierr); 342 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 343 CeedChk(ierr); 344 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 345 CeedChk(ierr); 346 ierr = CeedQFunctionFieldGetNumComponents(qfoutputfields[i], &ncomp); 347 CeedChk(ierr); 348 // Basis action 349 switch(emode) { 350 case CEED_EVAL_NONE: 351 break; // No action 352 case CEED_EVAL_INTERP: 353 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 354 CeedChk(ierr); 355 ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST, 356 CEED_USE_POINTER, 357 &impl->edata[i + numinputfields][e*elemsize*ncomp]); 358 ierr = CeedBasisApply(basis, blksize, CEED_TRANSPOSE, 359 CEED_EVAL_INTERP, impl->qvecsout[i], 360 impl->evecsout[i]); CeedChk(ierr); 361 break; 362 case CEED_EVAL_GRAD: 363 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 364 CeedChk(ierr); 365 ierr = CeedVectorSetArray(impl->evecsout[i], CEED_MEM_HOST, 366 CEED_USE_POINTER, 367 &impl->edata[i + numinputfields][e*elemsize*ncomp]); 368 ierr = CeedBasisApply(basis, blksize, CEED_TRANSPOSE, 369 CEED_EVAL_GRAD, impl->qvecsout[i], 370 impl->evecsout[i]); CeedChk(ierr); 371 break; 372 case CEED_EVAL_WEIGHT: { 373 Ceed ceed; 374 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 375 return CeedError(ceed, 1, 376 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 377 break; // Should not occur 378 } 379 case CEED_EVAL_DIV: 380 break; // Not implimented 381 case CEED_EVAL_CURL: 382 break; // Not implimented 383 } 384 } 385 } 386 387 // Zero lvecs 388 for (CeedInt i=0; i<numoutputfields; i++) { 389 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr); 390 if (vec == CEED_VECTOR_ACTIVE) 391 vec = outvec; 392 ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr); 393 } 394 395 // Output restriction 396 for (CeedInt i=0; i<numoutputfields; i++) { 397 // Restore evec 398 ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein], 399 &impl->edata[i + numinputfields]); CeedChk(ierr); 400 // Get output vector 401 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr); 402 // Active 403 if (vec == CEED_VECTOR_ACTIVE) 404 vec = outvec; 405 // Restrict 406 ierr = CeedOperatorFieldGetLMode(opoutputfields[i], &lmode); CeedChk(ierr); 407 ierr = CeedElemRestrictionApply(impl->blkrestr[i+impl->numein], CEED_TRANSPOSE, 408 lmode, impl->evecs[i+impl->numein], vec, 409 request); CeedChk(ierr); 410 411 } 412 413 // Restore input arrays 414 for (CeedInt i=0; i<numinputfields; i++) { 415 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 416 CeedChk(ierr); 417 if (emode == CEED_EVAL_WEIGHT) { // Skip 418 } else { 419 ierr = CeedVectorRestoreArrayRead(impl->evecs[i], 420 (const CeedScalar **) &impl->edata[i]); 421 CeedChk(ierr); 422 } 423 } 424 425 return 0; 426 } 427 428 int CeedOperatorCreate_Blocked(CeedOperator op) { 429 int ierr; 430 Ceed ceed; 431 ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr); 432 CeedOperator_Blocked *impl; 433 434 ierr = CeedCalloc(1, &impl); CeedChk(ierr); 435 ierr = CeedOperatorSetData(op, (void *)&impl); 436 437 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply", 438 CeedOperatorApply_Blocked); CeedChk(ierr); 439 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 440 CeedOperatorDestroy_Blocked); CeedChk(ierr); 441 return 0; 442 } 443