xref: /libCEED/backends/cuda-gen/ceed-cuda-gen-operator.c (revision 44554ea01e90fce366fc2a203c44be15754a38d6)
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   const CeedInt P1d = data->maxP1d;
119   const CeedInt thread1d = CeedIntMax(Q1d, P1d);
120   if (dim==1) {
121     const CeedInt elemsPerBlock = 32;
122     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
123                                            ? 1 : 0 );
124     CeedInt sharedMem = elemsPerBlock*thread1d*sizeof(CeedScalar);
125     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, 1,
126                                       elemsPerBlock, sharedMem, opargs);
127   } else if (dim==2) {
128     const CeedInt elemsPerBlock = thread1d<4? 16 : 2;
129     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
130                                            ? 1 : 0 );
131     CeedInt sharedMem = elemsPerBlock*thread1d*thread1d*sizeof(CeedScalar);
132     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, thread1d,
133                                       elemsPerBlock, sharedMem, opargs);
134   } else if (dim==3) {
135     const CeedInt elemsPerBlock = thread1d<6? 4 : (thread1d<8? 2 : 1);
136     CeedInt grid = nelem/elemsPerBlock + ( (nelem/elemsPerBlock*elemsPerBlock<nelem)
137                                            ? 1 : 0 );
138     CeedInt sharedMem = elemsPerBlock*thread1d*thread1d*sizeof(CeedScalar);
139     ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, thread1d, thread1d,
140                                       elemsPerBlock, sharedMem, opargs);
141   }
142   CeedChk(ierr);
143 
144   // Restore input arrays
145   for (CeedInt i = 0; i < numinputfields; i++) {
146     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
147     CeedChk(ierr);
148     if (emode == CEED_EVAL_WEIGHT) { // Skip
149     } else {
150       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
151       if (vec == CEED_VECTOR_ACTIVE) vec = invec;
152       ierr = CeedVectorRestoreArrayRead(vec, &data->fields.in[i]);
153       CeedChk(ierr);
154     }
155   }
156 
157   // Restore output arrays
158   for (CeedInt i = 0; i < numoutputfields; i++) {
159     ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
160     CeedChk(ierr);
161     if (emode == CEED_EVAL_WEIGHT) { // Skip
162     } else {
163       ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
164       if (vec == CEED_VECTOR_ACTIVE) vec = outvec;
165       // Check for multiple output modes
166       CeedInt index = -1;
167       for (CeedInt j = 0; j < i; j++) {
168         if (vec == outvecs[j]) {
169           index = j;
170           break;
171         }
172       }
173       if (index == -1) {
174         ierr = CeedVectorRestoreArray(vec, &data->fields.out[i]);
175         CeedChk(ierr);
176       }
177     }
178   }
179 
180   // Restore context data
181   if (ctx) {
182     ierr = CeedQFunctionContextRestoreData(ctx, &qf_data->d_c);
183     CeedChk(ierr);
184   }
185   return 0;
186 }
187 
188 //------------------------------------------------------------------------------
189 // Create FDM element inverse not supported
190 //------------------------------------------------------------------------------
191 static int CeedOperatorCreateFDMElementInverse_Cuda(CeedOperator op) {
192   // LCOV_EXCL_START
193   int ierr;
194   Ceed ceed;
195   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
196   return CeedError(ceed, 1, "Backend does not implement FDM inverse creation");
197   // LCOV_EXCL_STOP
198 }
199 
200 //------------------------------------------------------------------------------
201 // Create operator
202 //------------------------------------------------------------------------------
203 int CeedOperatorCreate_Cuda_gen(CeedOperator op) {
204   int ierr;
205   Ceed ceed;
206   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
207   CeedOperator_Cuda_gen *impl;
208 
209   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
210   ierr = CeedOperatorSetData(op, impl); CeedChk(ierr);
211 
212   ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
213                                 CeedOperatorCreateFDMElementInverse_Cuda);
214   CeedChk(ierr);
215   ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd",
216                                 CeedOperatorApplyAdd_Cuda_gen); CeedChk(ierr);
217   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
218                                 CeedOperatorDestroy_Cuda_gen); CeedChk(ierr);
219   return 0;
220 }
221 //------------------------------------------------------------------------------
222