xref: /libCEED/backends/cuda-gen/ceed-cuda-gen-operator.c (revision 68b2f3bb6a01fa010da0a56571daeae100b32294)
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-cuda-gen.h"
18 #include "ceed-cuda-gen-operator-build.h"
19 
20 //------------------------------------------------------------------------------
21 // Destroy operator
22 //------------------------------------------------------------------------------
23 static int CeedOperatorDestroy_Cuda_gen(CeedOperator op) {
24   int ierr;
25   CeedOperator_Cuda_gen *impl;
26   ierr = CeedOperatorGetData(op, &impl); CeedChk(ierr);
27   ierr = CeedFree(&impl); CeedChk(ierr);
28   return 0;
29 }
30 
31 //------------------------------------------------------------------------------
32 // Apply and add to output
33 //------------------------------------------------------------------------------
34 static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector invec,
35     CeedVector outvec, CeedRequest *request) {
36   int ierr;
37   Ceed ceed;
38   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
39   CeedOperator_Cuda_gen *data;
40   ierr = CeedOperatorGetData(op, &data); CeedChk(ierr);
41   CeedQFunction qf;
42   CeedQFunction_Cuda_gen *qf_data;
43   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
44   ierr = CeedQFunctionGetData(qf, &qf_data); CeedChk(ierr);
45   CeedInt nelem, numinputfields, numoutputfields;
46   ierr = CeedOperatorGetNumElements(op, &nelem); CeedChk(ierr);
47   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
48   CeedChk(ierr);
49   CeedOperatorField *opinputfields, *opoutputfields;
50   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
51   CeedChk(ierr);
52   CeedQFunctionField *qfinputfields, *qfoutputfields;
53   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
54   CeedChk(ierr);
55   CeedEvalMode emode;
56   CeedVector vec, outvecs[16] = {};
57 
58   //Creation of the operator
59   ierr = CeedCudaGenOperatorBuild(op); CeedChk(ierr);
60 
61   // Input vectors
62   for (CeedInt i = 0; i < numinputfields; i++) {
63     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
64     CeedChk(ierr);
65     if (emode == CEED_EVAL_WEIGHT) { // Skip
66       data->fields.in[i] = NULL;
67     } else {
68       // Get input vector
69       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
70       if (vec == CEED_VECTOR_ACTIVE) vec = invec;
71       ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->fields.in[i]);
72       CeedChk(ierr);
73     }
74   }
75 
76   // Output vectors
77   for (CeedInt i = 0; i < numoutputfields; i++) {
78     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
79     CeedChk(ierr);
80     if (emode == CEED_EVAL_WEIGHT) { // Skip
81       data->fields.out[i] = NULL;
82     } else {
83       // Get output vector
84       ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
85       if (vec == CEED_VECTOR_ACTIVE) vec = outvec;
86       outvecs[i] = vec;
87       // Check for multiple output modes
88       CeedInt index = -1;
89       for (CeedInt j = 0; j < i; j++) {
90         if (vec == outvecs[j]) {
91           index = j;
92           break;
93         }
94       }
95       if (index == -1) {
96         ierr = CeedVectorGetArray(vec, CEED_MEM_DEVICE, &data->fields.out[i]);
97         CeedChk(ierr);
98       } else {
99         data->fields.out[i] = data->fields.out[index];
100       }
101     }
102   }
103 
104   // Get context data
105   CeedQFunctionContext ctx;
106   ierr = CeedQFunctionGetInnerContext(qf, &ctx); CeedChk(ierr);
107   if (ctx) {
108     ierr = CeedQFunctionContextGetData(ctx, CEED_MEM_DEVICE, &qf_data->d_c);
109     CeedChk(ierr);
110   }
111 
112   // Apply operator
113   void *opargs[] = {(void *) &nelem, &qf_data->d_c, &data->indices,
114                     &data->fields, &data->B, &data->G, &data->W
115                    };
116   const CeedInt dim = data->dim;
117   const CeedInt Q1d = data->Q1d;
118   if (dim==1) {
119     const CeedInt elemsPerBlock = 32;
120     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
121                                            ? 1 : 0 );
122     CeedInt sharedMem = elemsPerBlock*Q1d*sizeof(CeedScalar);
123     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, Q1d, 1, elemsPerBlock,
124                                       sharedMem, opargs);
125   } else if (dim==2) {
126     const CeedInt elemsPerBlock = Q1d<4? 16 : 2;
127     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
128                                            ? 1 : 0 );
129     CeedInt sharedMem = elemsPerBlock*Q1d*Q1d*sizeof(CeedScalar);
130     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, Q1d, Q1d,
131                                       elemsPerBlock, sharedMem, opargs);
132   } else if (dim==3) {
133     const CeedInt elemsPerBlock = Q1d<6? 4 : (Q1d<8? 2 : 1);
134     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
135                                            ? 1 : 0 );
136     CeedInt sharedMem = elemsPerBlock*Q1d*Q1d*sizeof(CeedScalar);
137     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, Q1d, Q1d,
138                                       elemsPerBlock, sharedMem, opargs);
139   }
140   CeedChk(ierr);
141 
142   // Restore input arrays
143   for (CeedInt i = 0; i < numinputfields; i++) {
144     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
145     CeedChk(ierr);
146     if (emode == CEED_EVAL_WEIGHT) { // Skip
147     } else {
148       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
149       if (vec == CEED_VECTOR_ACTIVE) vec = invec;
150       ierr = CeedVectorRestoreArrayRead(vec, &data->fields.in[i]);
151       CeedChk(ierr);
152     }
153   }
154 
155   // Restore output arrays
156   for (CeedInt i = 0; i < numoutputfields; i++) {
157     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
158     CeedChk(ierr);
159     if (emode == CEED_EVAL_WEIGHT) { // Skip
160     } else {
161       ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
162       if (vec == CEED_VECTOR_ACTIVE) vec = outvec;
163       // Check for multiple output modes
164       CeedInt index = -1;
165       for (CeedInt j = 0; j < i; j++) {
166         if (vec == outvecs[j]) {
167           index = j;
168           break;
169         }
170       }
171       if (index == -1) {
172         ierr = CeedVectorRestoreArray(vec, &data->fields.out[i]);
173         CeedChk(ierr);
174       }
175     }
176   }
177 
178   // Restore context data
179   if (ctx) {
180     ierr = CeedQFunctionContextRestoreData(ctx, &qf_data->d_c);
181     CeedChk(ierr);
182   }
183   return 0;
184 }
185 
186 //------------------------------------------------------------------------------
187 // Create FDM element inverse not supported
188 //------------------------------------------------------------------------------
189 static int CeedOperatorCreateFDMElementInverse_Cuda(CeedOperator op) {
190   // LCOV_EXCL_START
191   int ierr;
192   Ceed ceed;
193   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
194   return CeedError(ceed, 1, "Backend does not implement FDM inverse creation");
195   // LCOV_EXCL_STOP
196 }
197 
198 //------------------------------------------------------------------------------
199 // Create operator
200 //------------------------------------------------------------------------------
201 int CeedOperatorCreate_Cuda_gen(CeedOperator op) {
202   int ierr;
203   Ceed ceed;
204   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
205   CeedOperator_Cuda_gen *impl;
206 
207   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
208   ierr = CeedOperatorSetData(op, impl); CeedChk(ierr);
209 
210   ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
211                                 CeedOperatorCreateFDMElementInverse_Cuda);
212   CeedChk(ierr);
213   ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd",
214                                 CeedOperatorApplyAdd_Cuda_gen); CeedChk(ierr);
215   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
216                                 CeedOperatorDestroy_Cuda_gen); CeedChk(ierr);
217   return 0;
218 }
219 //------------------------------------------------------------------------------
220