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