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 ierr = CeedVectorDestroy(&impl->etmp); CeedChk(ierr); 26 ierr = CeedVectorDestroy(&impl->qdata); CeedChk(ierr); 27 ierr = CeedFree(&op->data); CeedChk(ierr); 28 return 0; 29 } 30 31 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector qdata, 32 CeedVector ustate, 33 CeedVector residual, CeedRequest *request) { 34 CeedOperator_Ref *impl = op->data; 35 CeedVector etmp; 36 CeedInt Q; 37 const CeedInt nc = op->basis->ndof, dim = op->basis->dim; 38 CeedScalar *Eu; 39 char *qd; 40 int ierr; 41 CeedTransposeMode lmode = CEED_NOTRANSPOSE; 42 43 if (!impl->etmp) { 44 ierr = CeedVectorCreate(op->ceed, 45 nc * op->Erestrict->nelem * op->Erestrict->elemsize, 46 &impl->etmp); CeedChk(ierr); 47 // etmp is allocated when CeedVectorGetArray is called below 48 } 49 etmp = impl->etmp; 50 if (op->qf->inmode & ~CEED_EVAL_WEIGHT) { 51 ierr = CeedElemRestrictionApply(op->Erestrict, CEED_NOTRANSPOSE, 52 nc, lmode, ustate, etmp, 53 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 54 } 55 ierr = CeedBasisGetNumQuadraturePoints(op->basis, &Q); CeedChk(ierr); 56 ierr = CeedVectorGetArray(etmp, CEED_MEM_HOST, &Eu); CeedChk(ierr); 57 ierr = CeedVectorGetArray(qdata, CEED_MEM_HOST, (CeedScalar**)&qd); 58 CeedChk(ierr); 59 for (CeedInt e=0; e<op->Erestrict->nelem; e++) { 60 CeedScalar BEu[Q*nc*(dim+2)], BEv[Q*nc*(dim+2)], *out[5] = {0,0,0,0,0}; 61 const CeedScalar *in[5] = {0,0,0,0,0}; 62 // TODO: quadrature weights can be computed just once 63 ierr = CeedBasisApply(op->basis, CEED_NOTRANSPOSE, op->qf->inmode, 64 &Eu[e*op->Erestrict->elemsize*nc], BEu); 65 CeedChk(ierr); 66 CeedScalar *u_ptr = BEu, *v_ptr = BEv; 67 if (op->qf->inmode & CEED_EVAL_INTERP) { in[0] = u_ptr; u_ptr += Q*nc; } 68 if (op->qf->inmode & CEED_EVAL_GRAD) { in[1] = u_ptr; u_ptr += Q*nc*dim; } 69 if (op->qf->inmode & CEED_EVAL_WEIGHT) { in[4] = u_ptr; u_ptr += Q; } 70 if (op->qf->outmode & CEED_EVAL_INTERP) { out[0] = v_ptr; v_ptr += Q*nc; } 71 if (op->qf->outmode & CEED_EVAL_GRAD) { out[1] = v_ptr; v_ptr += Q*nc*dim; } 72 ierr = CeedQFunctionApply(op->qf, &qd[e*Q*op->qf->qdatasize], Q, in, out); 73 CeedChk(ierr); 74 ierr = CeedBasisApply(op->basis, CEED_TRANSPOSE, op->qf->outmode, BEv, 75 &Eu[e*op->Erestrict->elemsize*nc]); 76 CeedChk(ierr); 77 } 78 ierr = CeedVectorRestoreArray(etmp, &Eu); CeedChk(ierr); 79 ierr = CeedVectorRestoreArray(qdata, (CeedScalar**)&qd); CeedChk(ierr); 80 if (residual) { 81 CeedScalar *res; 82 ierr = CeedVectorGetArray(residual, CEED_MEM_HOST, &res); CeedChk(ierr); 83 for (int i = 0; i < residual->length; i++) 84 res[i] = (CeedScalar)0; 85 ierr = CeedElemRestrictionApply(op->Erestrict, CEED_TRANSPOSE, 86 nc, lmode, etmp, residual, 87 CEED_REQUEST_IMMEDIATE); CeedChk(ierr); 88 ierr = CeedVectorRestoreArray(residual, &res); CeedChk(ierr); 89 } 90 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) 91 *request = NULL; 92 return 0; 93 } 94 95 static int CeedOperatorGetQData_Ref(CeedOperator op, CeedVector *qdata) { 96 CeedOperator_Ref *impl = op->data; 97 int ierr; 98 99 if (!impl->qdata) { 100 CeedInt Q; 101 ierr = CeedBasisGetNumQuadraturePoints(op->basis, &Q); CeedChk(ierr); 102 ierr = CeedVectorCreate(op->ceed, 103 op->Erestrict->nelem * Q 104 * op->qf->qdatasize / sizeof(CeedScalar), 105 &impl->qdata); CeedChk(ierr); 106 } 107 *qdata = impl->qdata; 108 return 0; 109 } 110 111 int CeedOperatorCreate_Ref(CeedOperator op) { 112 CeedOperator_Ref *impl; 113 int ierr; 114 115 ierr = CeedCalloc(1, &impl); CeedChk(ierr); 116 op->data = impl; 117 op->Destroy = CeedOperatorDestroy_Ref; 118 op->Apply = CeedOperatorApply_Ref; 119 op->GetQData = CeedOperatorGetQData_Ref; 120 return 0; 121 } 122