1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #define CEED_DEBUG_COLOR 12
9
10 #include <ceed.h>
11 #include <ceed/backend.h>
12 #include <ceed/gen-tools.h>
13 #include <ceed/jit-tools.h>
14 #include <cuda_runtime.h>
15
16 #include <iostream>
17 #include <sstream>
18 #include <string>
19
20 #include "../cuda-ref/ceed-cuda-ref.h"
21 #include "../cuda-shared/ceed-cuda-shared.h"
22 #include "../cuda/ceed-cuda-common.h"
23 #include "../cuda/ceed-cuda-compile.h"
24 #include "ceed-cuda-gen.h"
25
26 struct FieldReuse_Cuda {
27 CeedInt index;
28 bool is_input;
29 CeedEvalMode eval_mode;
30 };
31
32 //------------------------------------------------------------------------------
33 // Determine type of operator
34 //------------------------------------------------------------------------------
CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed,CeedInt num_input_fields,CeedOperatorField * op_input_fields,CeedQFunctionField * qf_input_fields,CeedInt num_output_fields,CeedOperatorField * op_output_fields,CeedQFunctionField * qf_output_fields,CeedInt * max_P,CeedInt * max_P_1d,CeedInt * Q,CeedInt * Q_1d,CeedInt * max_dim,bool * is_all_tensor,bool * use_3d_slices)35 static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fields, CeedOperatorField *op_input_fields,
36 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields,
37 CeedQFunctionField *qf_output_fields, CeedInt *max_P, CeedInt *max_P_1d, CeedInt *Q, CeedInt *Q_1d,
38 CeedInt *max_dim, bool *is_all_tensor, bool *use_3d_slices) {
39 // Check if all are tensor
40 *is_all_tensor = true;
41 for (CeedInt i = 0; i < num_input_fields; i++) {
42 CeedBasis basis;
43
44 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
45 if (basis != CEED_BASIS_NONE) {
46 bool is_field_tensor;
47
48 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor));
49 *is_all_tensor = *is_all_tensor && is_field_tensor;
50 }
51 CeedCallBackend(CeedBasisDestroy(&basis));
52 }
53 for (CeedInt i = 0; i < num_output_fields; i++) {
54 CeedBasis basis;
55
56 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
57 if (basis != CEED_BASIS_NONE) {
58 bool is_field_tensor;
59
60 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor));
61 *is_all_tensor = *is_all_tensor && is_field_tensor;
62 }
63 CeedCallBackend(CeedBasisDestroy(&basis));
64 }
65
66 // Find max_P, max_P_1d, Q, and Q_1d
67 bool is_all_3d = true;
68
69 *max_P = 0;
70 *max_P_1d = 0;
71 *Q = 0;
72 *Q_1d = 0;
73 for (CeedInt i = 0; i < num_input_fields; i++) {
74 CeedBasis basis;
75
76 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
77 if (basis != CEED_BASIS_NONE) {
78 bool is_field_tensor;
79 CeedInt field_dim = 0, field_P = 0, field_P_1d = 0, field_Q = 0, field_Q_1d = 0;
80
81 // Check if 3D
82 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim));
83 is_all_3d = is_all_3d && (field_dim == 3);
84 *max_dim = CeedIntMax(*max_dim, field_dim);
85
86 // Collect P, P_1d, Q, and Q_1d
87 CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P));
88 *max_P = CeedIntMax(*max_P, field_P);
89 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor));
90 if (is_field_tensor) {
91 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d));
92 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d);
93 }
94 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q));
95 CeedCheck(*Q == 0 || field_Q == *Q, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible");
96 *Q = field_Q;
97 if (is_field_tensor) {
98 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d));
99 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible");
100 *Q_1d = field_Q_1d;
101 }
102 }
103 CeedCallBackend(CeedBasisDestroy(&basis));
104 }
105 for (CeedInt i = 0; i < num_output_fields; i++) {
106 CeedBasis basis;
107
108 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
109 if (basis != CEED_BASIS_NONE) {
110 bool is_field_tensor;
111 CeedInt field_dim = 0, field_P = 0, field_P_1d = 0, field_Q = 0, field_Q_1d = 0;
112
113 // Check if 3D
114 CeedCallBackend(CeedBasisGetDimension(basis, &field_dim));
115 is_all_3d = is_all_3d && (field_dim == 3);
116 *max_dim = CeedIntMax(*max_dim, field_dim);
117
118 // Collect P, P_1d, Q, and Q_1d
119 CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P));
120 *max_P = CeedIntMax(*max_P, field_P);
121 CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor));
122 if (is_field_tensor) {
123 CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d));
124 *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d);
125 }
126 CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q));
127 CeedCheck(*Q == 0 || field_Q == *Q, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible");
128 *Q = field_Q;
129 if (is_field_tensor) {
130 CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d));
131 CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible");
132 *Q_1d = field_Q_1d;
133 }
134 }
135 CeedCallBackend(CeedBasisDestroy(&basis));
136 }
137
138 // Only use 3D collocated gradient parallelization strategy when gradient is computed
139 *use_3d_slices = false;
140 if (is_all_3d && *is_all_tensor) {
141 bool was_grad_found = false;
142
143 for (CeedInt i = 0; i < num_input_fields; i++) {
144 CeedEvalMode eval_mode;
145
146 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
147 if (eval_mode == CEED_EVAL_GRAD) {
148 CeedBasis_Cuda_shared *basis_data;
149 CeedBasis basis;
150
151 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
152 CeedCallBackend(CeedBasisGetData(basis, &basis_data));
153 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true);
154 was_grad_found = true;
155 CeedCallBackend(CeedBasisDestroy(&basis));
156 }
157 }
158 for (CeedInt i = 0; i < num_output_fields; i++) {
159 CeedEvalMode eval_mode;
160
161 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
162 if (eval_mode == CEED_EVAL_GRAD) {
163 CeedBasis_Cuda_shared *basis_data;
164 CeedBasis basis;
165
166 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
167 CeedCallBackend(CeedBasisGetData(basis, &basis_data));
168 *use_3d_slices = basis_data->d_collo_grad_1d && (was_grad_found ? *use_3d_slices : true);
169 was_grad_found = true;
170 CeedCallBackend(CeedBasisDestroy(&basis));
171 }
172 }
173 }
174 return CEED_ERROR_SUCCESS;
175 }
176
177 //------------------------------------------------------------------------------
178 // Setup fields
179 //------------------------------------------------------------------------------
CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream & code,CeedOperator_Cuda_gen * data,Tab & tab,CeedInt i,CeedOperatorField op_field,CeedQFunctionField qf_field,FieldReuse_Cuda field_reuse,CeedInt max_dim,CeedInt Q,CeedInt Q_1d,bool is_input,bool is_all_tensor,bool is_at_points,bool use_3d_slices,bool skip_active_load)180 static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i,
181 CeedOperatorField op_field, CeedQFunctionField qf_field, FieldReuse_Cuda field_reuse,
182 CeedInt max_dim, CeedInt Q, CeedInt Q_1d, bool is_input, bool is_all_tensor, bool is_at_points,
183 bool use_3d_slices, bool skip_active_load) {
184 bool is_tensor = true, is_active = true;
185 CeedBasis basis;
186
187 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis));
188 if (basis != CEED_BASIS_NONE) CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
189 {
190 CeedVector vec;
191
192 CeedCallBackend(CeedOperatorFieldGetVector(op_field, &vec));
193 is_active = vec == CEED_VECTOR_ACTIVE;
194 CeedCallBackend(CeedVectorDestroy(&vec));
195 }
196
197 const char *field_name;
198 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i);
199 std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q";
200 std::string option_name = (is_input ? "inputs" : "outputs");
201 CeedEvalMode eval_mode = CEED_EVAL_NONE;
202 CeedInt elem_size = 0, num_comp = 0, dim = max_dim, P_1d = 0;
203 CeedElemRestriction elem_rstr;
204 CeedBasis_Cuda_shared *basis_data;
205
206 // Field reuse info
207 bool use_previous_field = field_reuse.index != -1;
208
209 CeedCallBackend(CeedOperatorFieldGetName(op_field, &field_name));
210 code << tab << "// -- " << (is_input ? "Input" : "Output") << " field " << i << ": " << field_name << "\n";
211
212 // Get field data
213 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr));
214 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) {
215 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
216 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
217 }
218 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
219 if (basis != CEED_BASIS_NONE) {
220 CeedCallBackend(CeedBasisGetData(basis, &basis_data));
221 CeedCallBackend(CeedBasisGetDimension(basis, &dim));
222 if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
223 else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d));
224 }
225 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode));
226
227 // Set field constants
228 code << tab << "const CeedInt dim" << var_suffix << " = " << dim << ";\n";
229 if (is_tensor && !is_all_tensor) {
230 CeedInt P = 0;
231
232 CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
233 code << tab << "const CeedInt P" << var_suffix << " = " << (basis == CEED_BASIS_NONE ? Q : P) << ";\n";
234 }
235 code << tab << "const CeedInt " << P_name << " = " << (basis == CEED_BASIS_NONE ? Q_1d : P_1d) << ";\n";
236 if (eval_mode != CEED_EVAL_WEIGHT) {
237 code << tab << "const CeedInt num_comp" << var_suffix << " = " << num_comp << ";\n";
238 }
239
240 // Load basis data
241 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
242 switch (eval_mode) {
243 case CEED_EVAL_NONE:
244 break;
245 case CEED_EVAL_INTERP:
246 if (is_at_points) {
247 // AtPoints
248 if (!basis_data->d_chebyshev_interp_1d) {
249 CeedSize interp_bytes;
250 CeedScalar *chebyshev_interp_1d;
251
252 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
253 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
254 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
255 CeedCallCuda(CeedBasisReturnCeed(basis), cudaMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes));
256 CeedCallCuda(CeedBasisReturnCeed(basis),
257 cudaMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
258 CeedCallBackend(CeedFree(&chebyshev_interp_1d));
259 }
260 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d;
261 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d;
262 } else {
263 // Standard quadrature
264 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d;
265 else data->B.outputs[i] = basis_data->d_interp_1d;
266 }
267 if (use_previous_field && !skip_active_load) {
268 std::string reuse_var = "s_B" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index));
269
270 code << tab << "CeedScalar *s_B" << var_suffix << " = " << reuse_var << ";\n";
271 } else {
272 bool is_collocated = false;
273
274 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated));
275 if ((is_active && skip_active_load) || (is_collocated && !is_at_points)) {
276 code << tab << "CeedScalar *s_B" << var_suffix << " = NULL;\n";
277 } else {
278 code << tab << "__shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n";
279 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n";
280 }
281 }
282 break;
283 case CEED_EVAL_GRAD:
284 if (is_at_points) {
285 // AtPoints
286 if (!basis_data->d_chebyshev_interp_1d) {
287 CeedSize interp_bytes;
288 CeedScalar *chebyshev_interp_1d;
289
290 interp_bytes = P_1d * Q_1d * sizeof(CeedScalar);
291 CeedCallBackend(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d));
292 CeedCallBackend(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d));
293 CeedCallCuda(CeedBasisReturnCeed(basis), cudaMalloc((void **)&basis_data->d_chebyshev_interp_1d, interp_bytes));
294 CeedCallCuda(CeedBasisReturnCeed(basis),
295 cudaMemcpy(basis_data->d_chebyshev_interp_1d, chebyshev_interp_1d, interp_bytes, cudaMemcpyHostToDevice));
296 CeedCallBackend(CeedFree(&chebyshev_interp_1d));
297 }
298 if (is_input) data->B.inputs[i] = basis_data->d_chebyshev_interp_1d;
299 else data->B.outputs[i] = basis_data->d_chebyshev_interp_1d;
300 } else {
301 // Standard quadrature
302 if (is_input) data->B.inputs[i] = basis_data->d_interp_1d;
303 else data->B.outputs[i] = basis_data->d_interp_1d;
304 }
305 if (is_tensor) {
306 if (use_previous_field && !skip_active_load) {
307 std::string reuse_var = "s_B" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index));
308
309 code << tab << "CeedScalar *s_B" << var_suffix << " = " << reuse_var << ";\n";
310 } else {
311 bool is_collocated = false;
312
313 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated));
314 if ((is_active && skip_active_load) || (is_collocated && !is_at_points)) {
315 code << tab << "CeedScalar *s_B" << var_suffix << " = NULL;\n";
316 } else {
317 code << tab << "__shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n";
318 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n";
319 }
320 }
321 }
322 if (is_at_points) break; // No G mat for AtPoints
323 if (use_3d_slices) {
324 if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d;
325 else data->G.outputs[i] = basis_data->d_collo_grad_1d;
326 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) {
327 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index));
328
329 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n";
330 } else if (is_active && skip_active_load) {
331 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n";
332 } else {
333 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n";
334 code << tab << "LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n";
335 }
336 } else {
337 bool has_collo_grad = basis_data->d_collo_grad_1d;
338
339 if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d;
340 else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d;
341 if (has_collo_grad) {
342 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) {
343 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index));
344
345 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n";
346 } else if (is_active && skip_active_load) {
347 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n";
348 } else {
349 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n";
350 code << tab << "LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n";
351 }
352 } else {
353 if (use_previous_field && field_reuse.eval_mode == CEED_EVAL_GRAD && !skip_active_load) {
354 std::string reuse_var = "s_G" + ((field_reuse.is_input ? "_in_" : "_out_") + std::to_string(field_reuse.index));
355
356 code << tab << "CeedScalar *s_G" << var_suffix << " = " << reuse_var << ";\n";
357 } else if (is_active && skip_active_load) {
358 code << tab << "CeedScalar *s_G" << var_suffix << " = NULL;\n";
359 } else {
360 code << tab << "__shared__ CeedScalar s_G" << var_suffix << "[" << P_name << "*" << Q_name << (is_tensor ? "" : "*dim")
361 << (is_tensor ? "" : var_suffix) << "];\n";
362 code << tab << "LoadMatrix<" << P_name << ", " << Q_name << (is_tensor ? "" : "*dim") << (is_tensor ? "" : var_suffix) << ">(data, G."
363 << option_name << "[" << i << "], s_G" << var_suffix << ");\n";
364 }
365 }
366 }
367 break;
368 case CEED_EVAL_WEIGHT:
369 break; // No action
370 // LCOV_EXCL_START
371 case CEED_EVAL_DIV:
372 case CEED_EVAL_CURL:
373 break; // TODO: Not implemented
374 // LCOV_EXCL_STOP
375 }
376 CeedCallBackend(CeedBasisDestroy(&basis));
377 return CEED_ERROR_SUCCESS;
378 }
379
380 //------------------------------------------------------------------------------
381 // Restriction
382 //------------------------------------------------------------------------------
CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream & code,CeedOperator_Cuda_gen * data,Tab & tab,CeedInt i,CeedInt field_input_buffer[],CeedOperatorField op_field,CeedQFunctionField qf_field,CeedInt max_dim,CeedInt Q_1d,bool is_input,bool is_all_tensor,bool is_at_points,bool use_3d_slices)383 static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i,
384 CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field,
385 CeedInt max_dim, CeedInt Q_1d, bool is_input, bool is_all_tensor, bool is_at_points,
386 bool use_3d_slices) {
387 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i);
388 std::string P_name = (is_all_tensor ? "P_1d" : "P") + var_suffix;
389 CeedEvalMode eval_mode = CEED_EVAL_NONE;
390 CeedInt elem_size = 0, num_comp = 0;
391 CeedSize l_size;
392 CeedRestrictionType rstr_type = CEED_RESTRICTION_STANDARD;
393 CeedElemRestriction_Cuda *rstr_data;
394 CeedElemRestriction elem_rstr;
395
396 // Get field data
397 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr));
398 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) {
399 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
400 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
401 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
402 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data));
403 }
404 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode));
405
406 // Restriction
407 if (is_input) {
408 // Input
409 if (field_input_buffer[i] != i) {
410 std::string buffer_name = "r_e_in_" + std::to_string(field_input_buffer[i]);
411
412 // Restriction was already done for previous input
413 code << tab << "CeedScalar *r_e" << var_suffix << " = " << buffer_name << ";\n";
414 } else if (eval_mode != CEED_EVAL_WEIGHT && !((eval_mode == CEED_EVAL_NONE) && use_3d_slices && is_at_points)) {
415 if (eval_mode == CEED_EVAL_NONE && rstr_type != CEED_RESTRICTION_POINTS) {
416 // No basis action, so r_e_in_* in also r_q_in_* and needs to be allocated
417 code << tab << "CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << P_name << "];\n";
418 } else if (rstr_type != CEED_RESTRICTION_POINTS) {
419 // Otherwise we're using the scratch space
420 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n";
421 }
422 switch (rstr_type) {
423 case CEED_RESTRICTION_STANDARD: {
424 CeedInt comp_stride;
425
426 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size));
427 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n";
428 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
429 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
430 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets;
431 code << tab << "ReadLVecStandard" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", "
432 << P_name << ">(data, l_size" << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix
433 << ");\n";
434 break;
435 }
436 case CEED_RESTRICTION_STRIDED: {
437 bool has_backend_strides;
438 CeedInt num_elem;
439
440 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides));
441 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem));
442 CeedInt strides[3] = {1, elem_size * num_elem, elem_size};
443
444 if (!has_backend_strides) {
445 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides));
446 }
447 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1]
448 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n";
449 code << tab << "ReadLVecStrided" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", " << P_name << ", strides"
450 << var_suffix << "_0, strides" << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, d" << var_suffix << ", r_e"
451 << var_suffix << ");\n";
452 break;
453 }
454 case CEED_RESTRICTION_POINTS: {
455 CeedInt comp_stride;
456
457 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
458 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
459 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets;
460 break;
461 }
462 // LCOV_EXCL_START
463 case CEED_RESTRICTION_ORIENTED:
464 case CEED_RESTRICTION_CURL_ORIENTED:
465 break; // TODO: Not implemented
466 // LCOV_EXCL_STOP
467 }
468 }
469 } else {
470 // Output
471 switch (rstr_type) {
472 case CEED_RESTRICTION_STANDARD: {
473 CeedInt comp_stride;
474
475 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size));
476 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n";
477 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
478 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
479 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets;
480 code << tab << "WriteLVecStandard" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", "
481 << P_name << ">(data, l_size" << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix
482 << ");\n";
483 break;
484 }
485 case CEED_RESTRICTION_STRIDED: {
486 bool has_backend_strides;
487 CeedInt num_elem;
488
489 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides));
490 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem));
491 CeedInt strides[3] = {1, elem_size * num_elem, elem_size};
492
493 if (!has_backend_strides) {
494 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides));
495 }
496 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1]
497 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n";
498 code << tab << "WriteLVecStrided" << (is_all_tensor ? max_dim : 1) << "d<num_comp" << var_suffix << ", " << P_name << ", strides"
499 << var_suffix << "_0, strides" << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, r_e" << var_suffix << ", d" << var_suffix
500 << ");\n";
501 break;
502 }
503 case CEED_RESTRICTION_POINTS:
504 data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets;
505 break;
506 // LCOV_EXCL_START
507 case CEED_RESTRICTION_ORIENTED:
508 case CEED_RESTRICTION_CURL_ORIENTED:
509 break; // TODO: Not implemented
510 // LCOV_EXCL_STOP
511 }
512 }
513 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
514 return CEED_ERROR_SUCCESS;
515 }
516
517 //------------------------------------------------------------------------------
518 // Basis
519 //------------------------------------------------------------------------------
CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream & code,CeedOperator_Cuda_gen * data,Tab & tab,CeedInt i,CeedOperatorField op_field,CeedQFunctionField qf_field,CeedInt max_dim,CeedInt Q_1d,bool is_input,bool is_all_tensor,bool is_at_points,bool use_3d_slices)520 static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt i,
521 CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt max_dim, CeedInt Q_1d,
522 bool is_input, bool is_all_tensor, bool is_at_points, bool use_3d_slices) {
523 bool is_tensor = true, is_collocated = true;
524 CeedBasis basis;
525 CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis));
526 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
527 CeedCallBackend(CeedBasisIsCollocated(basis, &is_collocated));
528
529 std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i);
530 std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q";
531 CeedEvalMode eval_mode = CEED_EVAL_NONE;
532 CeedInt dim = max_dim, elem_size = 0, num_comp = 0, P_1d = 0;
533 CeedElemRestriction elem_rstr;
534
535 // Get field data
536 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_field, &elem_rstr));
537 if (elem_rstr != CEED_ELEMRESTRICTION_NONE) {
538 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
539 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
540 }
541 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
542 if (basis != CEED_BASIS_NONE) {
543 CeedCallBackend(CeedBasisGetDimension(basis, &dim));
544 if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
545 else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d));
546 }
547 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode));
548
549 // Basis
550 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
551 if (is_input) {
552 switch (eval_mode) {
553 case CEED_EVAL_NONE:
554 if (!use_3d_slices && !is_at_points) {
555 code << tab << "CeedScalar *r_q" << var_suffix << " = r_e" << var_suffix << ";\n";
556 }
557 break;
558 case CEED_EVAL_INTERP:
559 if (is_at_points) {
560 std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d";
561
562 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n";
563 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix
564 << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n";
565 } else {
566 std::string function_name = is_tensor ? ((dim == 1 ? "Interp" : "InterpTensor") + std::string(is_collocated ? "CollocatedNodes" : "") +
567 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened"))
568 : "InterpNonTensor";
569 std::string op_t_1d_name = (is_all_tensor || !is_tensor) ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name);
570
571 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (dim >= 3) ? Q_name : "1") << "];\n";
572 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_e"
573 << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n";
574 }
575 break;
576 case CEED_EVAL_GRAD:
577 if (is_at_points) {
578 std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d";
579
580 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n";
581 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix
582 << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n";
583 } else if (use_3d_slices) {
584 std::string function_name =
585 (dim > 1 ? "InterpTensor" : "Interp") + std::string(is_collocated ? "CollocatedNodes" : "") + std::to_string(dim) + "d";
586
587 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n";
588 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_e" << var_suffix
589 << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n";
590 } else if (is_tensor) {
591 bool is_collocated_grad = dim == 3 && Q_1d >= P_1d;
592 std::string function_name =
593 (dim == 1 ? "Grad" : ("GradTensor" + std::string(is_collocated ? "CollocatedNodes" : (is_collocated_grad ? "Collocated" : "")))) +
594 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened");
595 std::string op_t_1d_name = is_all_tensor ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name);
596
597 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "*"
598 << (is_all_tensor && dim >= 3 ? Q_name : "1") << "];\n";
599 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_e"
600 << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n";
601 } else {
602 std::string function_name = "GradNonTensor";
603
604 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
605 code << tab << function_name << "<num_comp" << var_suffix << ", dim" << var_suffix << ", " << P_name << ", " << Q_name
606 << ", OP_T_1D>(data, r_e" << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n";
607 }
608 break;
609 case CEED_EVAL_WEIGHT: {
610 if (is_at_points) {
611 code << tab << "// Nothing to do AtPoints\n";
612 } else {
613 CeedBasis_Cuda_shared *basis_data;
614 std::string function_name = is_tensor
615 ? ((dim == 1 ? "Weight" : "WeightTensor") + std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened"))
616 : "WeightNonTensor";
617
618 code << tab << "CeedScalar r_q" << var_suffix << "[" << (is_all_tensor && (dim >= 3) ? Q_name : "1") << "];\n";
619 CeedCallBackend(CeedBasisGetData(basis, &basis_data));
620 data->W = basis_data->d_q_weight_1d;
621 code << tab << function_name << "<" << P_name << ", " << Q_name << ">(data, W, r_q" << var_suffix << ");\n";
622 }
623 break;
624 }
625 // LCOV_EXCL_START
626 case CEED_EVAL_DIV:
627 case CEED_EVAL_CURL:
628 break; // TODO: Not implemented
629 // LCOV_EXCL_STOP
630 }
631 } else {
632 switch (eval_mode) {
633 case CEED_EVAL_NONE:
634 code << tab << "CeedScalar *r_e" << var_suffix << " = r_q" << var_suffix << ";\n";
635 break; // No action
636 case CEED_EVAL_INTERP:
637 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n";
638 if (is_at_points) {
639 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d";
640
641 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_c" << var_suffix
642 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n";
643 } else {
644 std::string function_name =
645 is_tensor ? ((dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::string(is_collocated ? "CollocatedNodes" : "") +
646 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened"))
647 : "InterpTransposeNonTensor";
648 std::string op_t_1d_name = (is_all_tensor || !is_tensor) ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name);
649
650 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_q"
651 << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n";
652 }
653 break;
654 case CEED_EVAL_GRAD:
655 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n";
656 if (is_at_points) {
657 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d";
658
659 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_c" << var_suffix
660 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n";
661 } else if (use_3d_slices) {
662 std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::string(is_collocated ? "CollocatedNodes" : "") +
663 std::to_string(dim) + "d";
664
665 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", OP_T_1D>(data, r_q" << var_suffix
666 << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n";
667 } else if (is_tensor) {
668 bool is_collocated_grad = dim == 3 && Q_1d >= P_1d;
669 std::string function_name =
670 (dim == 1 ? "GradTranspose"
671 : ("GradTransposeTensor" + std::string(is_collocated ? "CollocatedNodes" : (is_collocated_grad ? "Collocated" : "")))) +
672 std::to_string(dim) + "d" + (is_all_tensor ? "" : "Flattened");
673 std::string op_t_1d_name = is_all_tensor ? "OP_T_1D" : (P_1d > Q_1d ? P_name : Q_name);
674
675 code << tab << function_name << "<num_comp" << var_suffix << ", " << P_name << ", " << Q_name << ", " << op_t_1d_name << ">(data, r_q"
676 << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n";
677 } else {
678 std::string function_name = "GradTransposeNonTensor";
679
680 code << tab << function_name << "<num_comp" << var_suffix << ", dim" << var_suffix << ", " << P_name << ", " << Q_name
681 << ", OP_T_1D>(data, r_q" << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n";
682 }
683 break;
684 // LCOV_EXCL_START
685 case CEED_EVAL_WEIGHT:
686 break; // Should not occur
687 case CEED_EVAL_DIV:
688 case CEED_EVAL_CURL:
689 break; // TODO: Not implemented
690 // LCOV_EXCL_STOP
691 }
692 }
693 CeedCallBackend(CeedBasisDestroy(&basis));
694 return CEED_ERROR_SUCCESS;
695 }
696
697 //------------------------------------------------------------------------------
698 // QFunction
699 //------------------------------------------------------------------------------
CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream & code,CeedOperator_Cuda_gen * data,Tab & tab,CeedInt max_dim,CeedInt max_num_points,CeedInt num_input_fields,CeedOperatorField * op_input_fields,CeedQFunctionField * qf_input_fields,CeedInt num_output_fields,CeedOperatorField * op_output_fields,CeedQFunctionField * qf_output_fields,std::string qfunction_name,CeedInt Q_1d,bool is_all_tensor,bool is_at_points,bool use_3d_slices,bool is_assemble)700 static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, Tab &tab, CeedInt max_dim,
701 CeedInt max_num_points, CeedInt num_input_fields, CeedOperatorField *op_input_fields,
702 CeedQFunctionField *qf_input_fields, CeedInt num_output_fields,
703 CeedOperatorField *op_output_fields, CeedQFunctionField *qf_output_fields,
704 std::string qfunction_name, CeedInt Q_1d, bool is_all_tensor, bool is_at_points,
705 bool use_3d_slices, bool is_assemble) {
706 std::string Q_name = is_all_tensor ? "Q_1d" : "Q";
707 CeedEvalMode eval_mode = CEED_EVAL_NONE;
708 CeedElemRestriction elem_rstr;
709
710 // Setup output arrays
711 code << "\n";
712 code << tab << "// -- Output field setup\n";
713 for (CeedInt i = 0; i < num_output_fields; i++) {
714 const char *field_name;
715 std::string var_suffix = "_out_" + std::to_string(i);
716
717 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
718 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
719 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
720 switch (eval_mode) {
721 case CEED_EVAL_NONE:
722 if (is_at_points) {
723 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n";
724 } else {
725 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (max_dim >= 3) ? Q_name : "1")
726 << "];\n";
727 }
728 break;
729 case CEED_EVAL_INTERP:
730 if (is_at_points) {
731 // Accumulator for point data
732 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "];\n";
733 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "; i++) r_c" << var_suffix
734 << "[i] = 0.0;\n";
735 } else {
736 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_all_tensor && (max_dim >= 3) ? Q_name : "1")
737 << "];\n";
738 }
739 break;
740 case CEED_EVAL_GRAD:
741 if (is_at_points) {
742 // Accumulator for point data
743 code << tab << "CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "];\n";
744 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (max_dim >= 3 ? Q_name : "1") << "; i++) r_c" << var_suffix
745 << "[i] = 0.0;\n";
746 } else if (use_3d_slices) {
747 // Accumulator for gradient slices
748 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n";
749 code << tab << "for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << Q_name << "; i++) r_q" << var_suffix << "[i] = 0.0;\n";
750 } else {
751 code << tab << "CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "*"
752 << (is_all_tensor && (max_dim >= 3) ? Q_name : "1") << "];\n";
753 }
754 break;
755 case CEED_EVAL_WEIGHT:
756 break;
757 // LCOV_EXCL_START
758 case CEED_EVAL_DIV:
759 case CEED_EVAL_CURL:
760 break; // TODO: Not implemented
761 // LCOV_EXCL_STOP
762 }
763 }
764
765 if (is_at_points) {
766 // We need to handle batches of points
767 code << "\n";
768 code << tab << "// Note: Using batches of points\n";
769 code << tab << "const CeedInt point_loop_bound = (blockDim.x*blockDim.y) * ceil((1.0*max_num_points) / (blockDim.x*blockDim.y));\n\n";
770 code << tab << "#pragma unroll\n";
771 code << tab << "for (CeedInt i = threadIdx.x + threadIdx.y*blockDim.x; i < point_loop_bound; i += blockDim.x*blockDim.y) {\n";
772 tab.push();
773 code << tab << "const CeedInt p = i % max_num_points;\n\n";
774
775 code << tab << "// -- Coordinates\n";
776 code << tab << "CeedScalar r_x[max_dim];\n";
777 code << tab << "ReadPoint<max_dim, coords_comp_stride, max_num_points>(data, elem, p, max_num_points, points.indices, points.coords, r_x);\n\n";
778
779 code << tab << "// -- Input fields\n";
780 for (CeedInt i = 0; i < num_input_fields; i++) {
781 const char *field_name;
782 std::string var_suffix = "_in_" + std::to_string(i);
783 std::string P_name = "P_1d" + var_suffix;
784
785 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name));
786 code << tab << "// ---- Input field " << i << ": " << field_name << "\n";
787 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
788 // Basis action
789 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
790 switch (eval_mode) {
791 case CEED_EVAL_NONE:
792 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
793 code << tab << "ReadPoint<num_comp" << var_suffix << ", comp_stride" << var_suffix
794 << ", max_num_points>(data, elem, p, max_num_points, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n";
795 break;
796 case CEED_EVAL_INTERP:
797 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
798 code << tab << "InterpAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
799 << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n";
800 break;
801 case CEED_EVAL_GRAD:
802 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
803 code << tab << "GradAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
804 << ">(data, i, r_c" << var_suffix << ", r_x, r_s" << var_suffix << ");\n";
805 break;
806 case CEED_EVAL_WEIGHT:
807 code << tab << "CeedScalar r_s" << var_suffix << "[1];\n";
808 code << tab << "r_s" << var_suffix << "[0] = 1.0;\n";
809 break;
810 // LCOV_EXCL_START
811 case CEED_EVAL_DIV:
812 case CEED_EVAL_CURL:
813 break; // TODO: Not implemented
814 // LCOV_EXCL_STOP
815 }
816 }
817 code << "\n";
818 code << tab << "// -- Output fields\n";
819 for (CeedInt i = 0; i < num_output_fields; i++) {
820 const char *field_name;
821 std::string var_suffix = "_out_" + std::to_string(i);
822
823 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
824 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
825 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
826 // Basis action
827 switch (eval_mode) {
828 case CEED_EVAL_NONE:
829 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
830 break;
831 case CEED_EVAL_INTERP:
832 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
833 break;
834 case CEED_EVAL_GRAD:
835 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
836 break;
837 // LCOV_EXCL_START
838 case CEED_EVAL_WEIGHT:
839 break; // Should not occur
840 case CEED_EVAL_DIV:
841 case CEED_EVAL_CURL:
842 break; // TODO: Not implemented
843 // LCOV_EXCL_STOP
844 }
845 }
846
847 } else if (use_3d_slices) {
848 // We treat quadrature points per slice in 3d to save registers
849 code << "\n";
850 code << tab << "// Note: Using planes of 3D elements\n";
851 code << tab << "#pragma unroll\n";
852 code << tab << "for (CeedInt q = 0; q < " << Q_name << "; q++) {\n";
853 tab.push();
854 code << tab << "// -- Input fields\n";
855 for (CeedInt i = 0; i < num_input_fields; i++) {
856 const char *field_name;
857 std::string var_suffix = "_in_" + std::to_string(i);
858
859 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name));
860 code << tab << "// ---- Input field " << i << ": " << field_name << "\n";
861 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
862 // Basis action
863 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
864 switch (eval_mode) {
865 case CEED_EVAL_NONE:
866 bool is_strided;
867
868 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
869
870 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
871 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided));
872 if (is_strided) {
873 bool has_backend_strides;
874 CeedInt num_elem, elem_size;
875
876 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
877 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &has_backend_strides));
878 CeedCallBackend(CeedElemRestrictionGetNumElements(elem_rstr, &num_elem));
879 CeedInt strides[3] = {1, elem_size * num_elem, elem_size};
880
881 if (!has_backend_strides) {
882 CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides));
883 }
884 code << tab << "const CeedInt strides" << var_suffix << "_0 = " << strides[0] << ", strides" << var_suffix << "_1 = " << strides[1]
885 << ", strides" << var_suffix << "_2 = " << strides[2] << ";\n";
886 code << tab << "ReadEVecSliceStrided3d<num_comp" << var_suffix << ", " << Q_name << ", strides" << var_suffix << "_0, strides"
887 << var_suffix << "_1, strides" << var_suffix << "_2>(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n";
888 } else {
889 CeedSize l_size = 0;
890 CeedInt comp_stride;
891 CeedElemRestriction_Cuda *rstr_data;
892
893 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size));
894 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n";
895 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
896 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
897 CeedCallBackend(CeedElemRestrictionGetData(elem_rstr, &rstr_data));
898 data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets;
899 code << tab << "ReadEVecSliceStandard3d<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", " << Q_name << ">(data, l_size"
900 << var_suffix << ", elem, q, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n";
901 }
902 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
903 break;
904 case CEED_EVAL_INTERP:
905 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
906 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) {\n";
907 tab.push();
908 code << tab << "r_s" << var_suffix << "[j] = r_q" << var_suffix << "[q + j*" << Q_name << "];\n";
909 tab.pop();
910 code << tab << "}\n";
911 break;
912 case CEED_EVAL_GRAD:
913 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
914 code << tab << "GradColloSlice3d<num_comp" << var_suffix << ", " << Q_name << ", OP_T_1D>(data, q, r_q" << var_suffix << ", s_G"
915 << var_suffix << ", r_s" << var_suffix << ");\n";
916 break;
917 case CEED_EVAL_WEIGHT:
918 code << tab << "CeedScalar r_s" << var_suffix << "[1];\n";
919 code << tab << "r_s" << var_suffix << "[0] = r_q" << var_suffix << "[q];\n";
920 break;
921 // LCOV_EXCL_START
922 case CEED_EVAL_DIV:
923 case CEED_EVAL_CURL:
924 break; // TODO: Not implemented
925 // LCOV_EXCL_STOP
926 }
927 }
928 code << "\n";
929 code << tab << "// -- Output fields\n";
930 for (CeedInt i = 0; i < num_output_fields; i++) {
931 const char *field_name;
932 std::string var_suffix = "_out_" + std::to_string(i);
933
934 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
935 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
936 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
937 // Basis action
938 switch (eval_mode) {
939 case CEED_EVAL_NONE:
940 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
941 break;
942 case CEED_EVAL_INTERP:
943 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n";
944 break;
945 case CEED_EVAL_GRAD:
946 code << tab << "CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "*dim" << var_suffix << "];\n";
947 break;
948 // LCOV_EXCL_START
949 case CEED_EVAL_WEIGHT:
950 break; // Should not occur
951 case CEED_EVAL_DIV:
952 case CEED_EVAL_CURL:
953 break; // TODO: Not implemented
954 // LCOV_EXCL_STOP
955 }
956 }
957 } else {
958 code << "\n";
959 code << tab << "// Note: Using full elements\n";
960 code << tab << "{\n";
961 tab.push();
962 code << tab << "// -- Input fields\n";
963 for (CeedInt i = 0; i < num_input_fields; i++) {
964 const char *field_name;
965
966 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name));
967 code << tab << "// ---- Input field " << i << ": " << field_name << "\n";
968 code << tab << "CeedScalar *r_s_in_" << i << " = r_q_in_" << i << ";\n";
969 }
970 code << tab << "// -- Output fields\n";
971 for (CeedInt i = 0; i < num_output_fields; i++) {
972 const char *field_name;
973
974 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
975 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
976 code << tab << "CeedScalar *r_s_out_" << i << " = r_q_out_" << i << ";\n";
977 }
978 }
979
980 // Input and output buffers
981 code << "\n";
982 code << tab << "// -- QFunction inputs and outputs\n";
983 code << tab << "// ---- Inputs\n";
984 code << tab << "CeedScalar *inputs[" << CeedIntMax(num_input_fields, 1) << "];\n";
985 for (CeedInt i = 0; i < num_input_fields; i++) {
986 const char *field_name;
987
988 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[i], &field_name));
989 code << tab << "// ------ Input field " << i << ": " << field_name << "\n";
990 code << tab << "inputs[" << i << "] = r_s_in_" << i << ";\n";
991 }
992 code << tab << "// ---- Outputs\n";
993 code << tab << "CeedScalar *outputs[" << CeedIntMax(num_output_fields, 1) << "];\n";
994 for (CeedInt i = 0; i < num_output_fields; i++) {
995 const char *field_name;
996
997 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
998 code << tab << "// ------ Output field " << i << ": " << field_name << "\n";
999 code << tab << "outputs[" << i << "] = r_s_out_" << i << ";\n";
1000 }
1001
1002 // Apply QFunction
1003 code << "\n";
1004 code << tab << "// -- Apply QFunction\n";
1005 code << tab << "" << qfunction_name << "(ctx, ";
1006 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) {
1007 code << "1";
1008 } else {
1009 code << Q_name;
1010 }
1011 code << ", inputs, outputs);\n";
1012
1013 if (is_at_points) {
1014 // Map back to coefficients
1015 code << "\n";
1016 code << tab << "// -- Output fields\n";
1017 for (CeedInt i = 0; i < num_output_fields; i++) {
1018 const char *field_name;
1019 std::string var_suffix = "_out_" + std::to_string(i);
1020 std::string P_name = "P_1d" + var_suffix;
1021
1022 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
1023 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
1024 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
1025 // Basis action
1026 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
1027 switch (eval_mode) {
1028 case CEED_EVAL_NONE: {
1029 CeedInt comp_stride;
1030 CeedElemRestriction elem_rstr;
1031
1032 if (is_assemble) break;
1033 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
1034 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
1035 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1036 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
1037 code << tab << "WritePoint<num_comp" << var_suffix << ", comp_stride" << var_suffix
1038 << ", max_num_points>(data, elem, i, points.num_per_elem[elem], indices.outputs[" << i << "]"
1039 << ", r_s" << var_suffix << ", d" << var_suffix << ");\n";
1040 break;
1041 }
1042 case CEED_EVAL_INTERP:
1043 code << tab << "if (i >= points.num_per_elem[elem]) {\n";
1044 tab.push();
1045 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n";
1046 tab.pop();
1047 code << tab << "}\n";
1048 code << tab << "InterpTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
1049 << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
1050 break;
1051 case CEED_EVAL_GRAD:
1052 code << tab << "if (i >= points.num_per_elem[elem]) {\n";
1053 tab.push();
1054 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << "*dim" << var_suffix << "; j++) r_s" << var_suffix << "[j] = 0.0;\n";
1055 tab.pop();
1056 code << tab << "}\n";
1057 code << tab << "GradTransposeAtPoints" << max_dim << "d<num_comp" << var_suffix << ", max_num_points, " << P_name << ", " << Q_name
1058 << ">(data, i, r_s" << var_suffix << ", r_x, r_c" << var_suffix << ");\n";
1059 break;
1060 // LCOV_EXCL_START
1061 case CEED_EVAL_WEIGHT:
1062 break; // Should not occur
1063 case CEED_EVAL_DIV:
1064 case CEED_EVAL_CURL:
1065 break; // TODO: Not implemented
1066 // LCOV_EXCL_STOP
1067 }
1068 }
1069 } else if (use_3d_slices) {
1070 // Copy or apply transpose grad, if needed
1071 code << "\n";
1072 code << tab << "// -- Output fields\n";
1073 for (CeedInt i = 0; i < num_output_fields; i++) {
1074 const char *field_name;
1075 std::string var_suffix = "_out_" + std::to_string(i);
1076 std::string P_name = "P_1d" + var_suffix;
1077
1078 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
1079 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
1080 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
1081 // Basis action
1082 code << tab << "// EvalMode: " << CeedEvalModes[eval_mode] << "\n";
1083 switch (eval_mode) {
1084 case CEED_EVAL_NONE:
1085 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n";
1086 tab.push();
1087 code << tab << "r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n";
1088 tab.pop();
1089 code << tab << "}\n";
1090 break;
1091 case CEED_EVAL_INTERP:
1092 code << tab << "for (CeedInt j = 0; j < num_comp" << var_suffix << " ; j++) {\n";
1093 tab.push();
1094 code << tab << "r_q" << var_suffix << "[q + j*" << Q_name << "] = r_s" << var_suffix << "[j];\n";
1095 tab.pop();
1096 code << tab << "}\n";
1097 break;
1098 case CEED_EVAL_GRAD:
1099 code << tab << "GradColloSliceTranspose3d<num_comp" << var_suffix << ", " << Q_name << ", OP_T_1D>(data, q, r_s" << var_suffix << ", s_G"
1100 << var_suffix << ", r_q" << var_suffix << ");\n";
1101 break;
1102 // LCOV_EXCL_START
1103 case CEED_EVAL_WEIGHT:
1104 break; // Should not occur
1105 case CEED_EVAL_DIV:
1106 case CEED_EVAL_CURL:
1107 break; // TODO: Not implemented
1108 // LCOV_EXCL_STOP
1109 }
1110 }
1111 }
1112 tab.pop();
1113 code << tab << "}\n";
1114 return CEED_ERROR_SUCCESS;
1115 }
1116
1117 //------------------------------------------------------------------------------
1118 // Build single operator kernel
1119 //------------------------------------------------------------------------------
CeedOperatorBuildKernel_Cuda_gen(CeedOperator op,bool * is_good_build)1120 extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op, bool *is_good_build) {
1121 bool is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false;
1122 Ceed ceed;
1123 CeedInt Q = 0, Q_1d = 0, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0;
1124 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1125 CeedQFunction_Cuda_gen *qf_data;
1126 CeedQFunction qf;
1127 CeedOperatorField *op_input_fields, *op_output_fields;
1128 CeedOperator_Cuda_gen *data;
1129 std::ostringstream code;
1130 Tab tab;
1131
1132 CeedCallBackend(CeedOperatorGetData(op, &data));
1133 {
1134 bool is_setup_done;
1135
1136 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done));
1137 if (is_setup_done) {
1138 *is_good_build = !data->use_fallback;
1139 return CEED_ERROR_SUCCESS;
1140 }
1141 }
1142
1143 // Check field compatibility
1144 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1145 {
1146 bool has_shared_bases = true;
1147
1148 for (CeedInt i = 0; i < num_input_fields; i++) {
1149 CeedBasis basis;
1150
1151 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
1152 if (basis != CEED_BASIS_NONE) {
1153 bool is_tensor = true;
1154 const char *resource;
1155 char *resource_root;
1156 Ceed basis_ceed;
1157
1158 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
1159 is_all_tensor = is_all_tensor && is_tensor;
1160 is_all_nontensor = is_all_nontensor && !is_tensor;
1161 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed));
1162 CeedCallBackend(CeedGetResource(basis_ceed, &resource));
1163 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root));
1164 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared");
1165 CeedCallBackend(CeedFree(&resource_root));
1166 CeedCallBackend(CeedDestroy(&basis_ceed));
1167 }
1168 CeedCallBackend(CeedBasisDestroy(&basis));
1169 }
1170
1171 for (CeedInt i = 0; i < num_output_fields; i++) {
1172 CeedBasis basis;
1173
1174 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
1175 if (basis != CEED_BASIS_NONE) {
1176 bool is_tensor = true;
1177 const char *resource;
1178 char *resource_root;
1179 Ceed basis_ceed;
1180
1181 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
1182 is_all_tensor = is_all_tensor && is_tensor;
1183 is_all_nontensor = is_all_nontensor && !is_tensor;
1184
1185 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed));
1186 CeedCallBackend(CeedGetResource(basis_ceed, &resource));
1187 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root));
1188 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared");
1189 CeedCallBackend(CeedFree(&resource_root));
1190 CeedCallBackend(CeedDestroy(&basis_ceed));
1191 }
1192 CeedCallBackend(CeedBasisDestroy(&basis));
1193 }
1194 // -- Fallback to ref if not all bases are shared
1195 if (!has_shared_bases) {
1196 *is_good_build = false;
1197 return CEED_ERROR_SUCCESS;
1198 }
1199 }
1200 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1201 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1202 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data));
1203 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1204
1205 // Get operator data
1206 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points));
1207 {
1208 CeedInt max_P = 0, max_P_1d = 0;
1209
1210 CeedCallBackend(CeedOperatorBuildKernelData_Cuda_gen(ceed, num_input_fields, op_input_fields, qf_input_fields, num_output_fields,
1211 op_output_fields, qf_output_fields, &max_P, &max_P_1d, &Q, &Q_1d, &max_dim, &is_all_tensor,
1212 &use_3d_slices));
1213 data->max_P_1d = is_all_tensor ? max_P_1d : max_P;
1214 }
1215 if (is_at_points) {
1216 CeedInt coords_dim = 0;
1217 CeedElemRestriction_Cuda *rstr_data;
1218 CeedElemRestriction rstr_points = NULL;
1219
1220 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
1221 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
1222 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride));
1223 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &coords_dim));
1224 CeedCallBackend(CeedElemRestrictionGetData(rstr_points, &rstr_data));
1225 data->points.indices = (CeedInt *)rstr_data->d_offsets;
1226 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1227 if (max_dim == 0) max_dim = coords_dim;
1228 if (Q_1d == 0) max_num_points = ceil(pow(max_num_points, 1.0 / max_dim));
1229 }
1230 if (max_dim == 0) max_dim = 1;
1231 data->dim = max_dim;
1232 if (is_at_points) use_3d_slices = false;
1233 if (Q_1d == 0) {
1234 if (is_at_points) Q_1d = max_num_points;
1235 else CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q_1d));
1236 }
1237 if (Q == 0) Q = Q_1d;
1238 data->Q = Q;
1239 data->Q_1d = Q_1d;
1240
1241 // Check for restriction only identity operator
1242 {
1243 bool is_identity_qf;
1244
1245 CeedCallBackend(CeedQFunctionIsIdentity(qf, &is_identity_qf));
1246 if (is_identity_qf) {
1247 CeedEvalMode eval_mode_in, eval_mode_out;
1248
1249 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[0], &eval_mode_in));
1250 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[0], &eval_mode_out));
1251 CeedCheck(eval_mode_in != CEED_EVAL_NONE || eval_mode_out != CEED_EVAL_NONE, ceed, CEED_ERROR_BACKEND,
1252 "Backend does not implement restriction only identity operators");
1253 }
1254 }
1255
1256 // Add atomicAdd function for old NVidia architectures
1257 {
1258 Ceed_Cuda *ceed_data;
1259 struct cudaDeviceProp prop;
1260
1261 CeedCallBackend(CeedGetData(ceed, &ceed_data));
1262 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id));
1263 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) {
1264 code << tab << "// AtomicAdd fallback source\n";
1265 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n";
1266 }
1267 }
1268
1269 // Load basis source files
1270 if (!is_all_nontensor) {
1271 code << tab << "// Tensor basis source\n";
1272 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n";
1273 }
1274 if (!is_all_tensor) {
1275 code << tab << "// Non-tensor basis source\n";
1276 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-nontensor-templates.h>\n\n";
1277 }
1278 if (!is_all_tensor && !is_all_nontensor) {
1279 code << "// Tensor basis source\n";
1280 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-flattened-templates.h>\n\n";
1281 }
1282 if (is_at_points) {
1283 code << "// AtPoints basis source\n";
1284 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h>\n\n";
1285 }
1286 code << "// CodeGen operator source\n";
1287 code << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n";
1288
1289 // Get QFunction name
1290 std::string qfunction_name(qf_data->qfunction_name);
1291 std::string operator_name;
1292
1293 operator_name = "CeedKernelCudaGenOperator_" + qfunction_name;
1294
1295 // Define CEED_Q_VLA
1296 code << "\n" << tab << "#undef CEED_Q_VLA\n";
1297 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) {
1298 code << tab << "#define CEED_Q_VLA 1\n\n";
1299 } else {
1300 code << tab << "#define CEED_Q_VLA " << Q_1d << "\n\n";
1301 }
1302
1303 // Add user QFunction source
1304 {
1305 const char *source_path;
1306
1307 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path));
1308 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file");
1309
1310 code << tab << "// User QFunction source\n";
1311 code << tab << "#include \"" << source_path << "\"\n\n";
1312 }
1313
1314 // Setup
1315 code << "\n" << tab << "// -----------------------------------------------------------------------------\n";
1316 code << tab << "// Operator Kernel\n";
1317 code << tab << "// \n";
1318 code << tab << "// d_[in,out]_i: CeedVector device array\n";
1319 code << tab << "// r_[in,out]_e_i: Element vector register\n";
1320 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n";
1321 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n";
1322 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n";
1323 code << tab << "// \n";
1324 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n";
1325 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n";
1326 code << tab << "// -----------------------------------------------------------------------------\n";
1327 code << tab << "extern \"C\" __global__ void " << operator_name
1328 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda "
1329 "points) {\n";
1330 tab.push();
1331
1332 // Scratch buffers
1333 for (CeedInt i = 0; i < num_input_fields; i++) {
1334 CeedEvalMode eval_mode;
1335
1336 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1337 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT
1338 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n";
1339 }
1340 }
1341 for (CeedInt i = 0; i < num_output_fields; i++) {
1342 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n";
1343 }
1344
1345 code << tab << "const CeedInt max_dim = " << max_dim << ";\n";
1346 if (!is_all_tensor) {
1347 code << tab << "const CeedInt Q = " << Q << ";\n";
1348 }
1349 if (!is_all_nontensor) {
1350 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n";
1351 }
1352 if (is_at_points) {
1353 code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n";
1354 code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n";
1355 }
1356
1357 // Shared data
1358 code << tab << "extern __shared__ CeedScalar slice[];\n";
1359 code << tab << "SharedData_Cuda data;\n";
1360 code << tab << "data.t_id_x = threadIdx.x;\n";
1361 code << tab << "data.t_id_y = threadIdx.y;\n";
1362 code << tab << "data.t_id_z = threadIdx.z;\n";
1363 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n";
1364 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n";
1365
1366 // -- Determine input mat reuse
1367 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX];
1368
1369 for (CeedInt i = 0; i < num_input_fields; i++) {
1370 input_matrix_reuse[i].index = -1;
1371 }
1372 for (CeedInt i = 0; i < num_input_fields; i++) {
1373 bool is_tensor = true;
1374 CeedEvalMode eval_mode_i;
1375 CeedBasis basis_i;
1376
1377 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i));
1378 if (eval_mode_i == CEED_EVAL_WEIGHT) continue;
1379 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i));
1380 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor));
1381 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) {
1382 CeedEvalMode eval_mode_j;
1383 CeedBasis basis_j;
1384
1385 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
1386 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1387 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
1388 if (basis_i == basis_j) {
1389 if (is_tensor) {
1390 input_matrix_reuse[i].index = j;
1391 input_matrix_reuse[i].is_input = true;
1392 input_matrix_reuse[i].eval_mode = eval_mode_j;
1393 } else {
1394 // For non-tensor can only re-use with the same eval mode
1395 if (eval_mode_i == eval_mode_j) {
1396 input_matrix_reuse[i].index = j;
1397 input_matrix_reuse[i].is_input = true;
1398 input_matrix_reuse[i].eval_mode = eval_mode_j;
1399 }
1400 }
1401 }
1402 CeedCallBackend(CeedBasisDestroy(&basis_j));
1403 }
1404 CeedCallBackend(CeedBasisDestroy(&basis_i));
1405 }
1406
1407 // -- Determine output mat reuse
1408 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX];
1409
1410 for (CeedInt i = 0; i < num_output_fields; i++) {
1411 output_matrix_reuse[i].index = -1;
1412 }
1413 for (CeedInt i = 0; i < num_output_fields; i++) {
1414 bool is_tensor = true;
1415 CeedEvalMode eval_mode_i;
1416 CeedBasis basis_i;
1417
1418 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i));
1419 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i));
1420 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor));
1421 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) {
1422 CeedEvalMode eval_mode_j;
1423 CeedBasis basis_j;
1424
1425 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
1426 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1427 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
1428 if (basis_i == basis_j) {
1429 if (is_tensor) {
1430 output_matrix_reuse[i].index = j;
1431 output_matrix_reuse[i].is_input = true;
1432 output_matrix_reuse[i].eval_mode = eval_mode_j;
1433 } else {
1434 // For non-tensor can only re-use with the same eval mode
1435 if (eval_mode_i == eval_mode_j) {
1436 output_matrix_reuse[i].index = j;
1437 output_matrix_reuse[i].is_input = true;
1438 output_matrix_reuse[i].eval_mode = eval_mode_j;
1439 }
1440 }
1441 }
1442 CeedCallBackend(CeedBasisDestroy(&basis_j));
1443 }
1444 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) {
1445 CeedEvalMode eval_mode_j;
1446 CeedBasis basis_j;
1447
1448 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j));
1449 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1450 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j));
1451 if (basis_i == basis_j) {
1452 if (is_tensor) {
1453 output_matrix_reuse[i].index = j;
1454 output_matrix_reuse[i].is_input = false;
1455 output_matrix_reuse[i].eval_mode = eval_mode_j;
1456 } else {
1457 // For non-tensor can only re-use with the same eval mode
1458 if (eval_mode_i == eval_mode_j) {
1459 output_matrix_reuse[i].index = j;
1460 output_matrix_reuse[i].is_input = false;
1461 output_matrix_reuse[i].eval_mode = eval_mode_j;
1462 }
1463 }
1464 }
1465 CeedCallBackend(CeedBasisDestroy(&basis_j));
1466 }
1467 CeedCallBackend(CeedBasisDestroy(&basis_i));
1468 }
1469
1470 // Initialize constants, and matrices B and G
1471 code << "\n" << tab << "// Input field constants and basis data\n";
1472 for (CeedInt i = 0; i < num_input_fields; i++) {
1473 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i],
1474 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, false));
1475 }
1476 code << "\n" << tab << "// Output field constants and basis data\n";
1477 for (CeedInt i = 0; i < num_output_fields; i++) {
1478 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i],
1479 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, false));
1480 }
1481
1482 // Loop over all elements
1483 code << "\n" << tab << "// Element loop\n";
1484 code << tab << "__syncthreads();\n";
1485 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n";
1486 tab.push();
1487
1488 // -- Compute minimum buffer space needed
1489 CeedInt max_rstr_buffer_size = 1;
1490
1491 for (CeedInt i = 0; i < num_input_fields; i++) {
1492 CeedEvalMode eval_mode;
1493
1494 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1495 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) {
1496 CeedInt num_comp;
1497 CeedElemRestriction elem_rstr;
1498
1499 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1500 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1501 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
1502 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1503 }
1504 }
1505 for (CeedInt i = 0; i < num_output_fields; i++) {
1506 CeedEvalMode eval_mode;
1507
1508 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
1509 if (eval_mode != CEED_EVAL_NONE) {
1510 CeedInt num_comp;
1511 CeedElemRestriction elem_rstr;
1512
1513 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
1514 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1515 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
1516 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1517 }
1518 }
1519 code << tab << "// Scratch restriction buffer space\n";
1520 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n";
1521
1522 // -- Determine best input field processing order
1523 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX];
1524
1525 for (CeedInt i = 0; i < num_input_fields; i++) {
1526 field_rstr_in_buffer[i] = -1;
1527 input_field_order[i] = -1;
1528 }
1529 {
1530 bool is_ordered[CEED_FIELD_MAX];
1531 CeedInt curr_index = 0;
1532
1533 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false;
1534 for (CeedInt i = 0; i < num_input_fields; i++) {
1535 CeedVector vec_i;
1536 CeedElemRestriction rstr_i;
1537
1538 if (is_ordered[i]) continue;
1539 field_rstr_in_buffer[i] = i;
1540 is_ordered[i] = true;
1541 input_field_order[curr_index] = i;
1542 curr_index++;
1543 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i));
1544 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT
1545 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i));
1546 for (CeedInt j = i + 1; j < num_input_fields; j++) {
1547 CeedVector vec_j;
1548 CeedElemRestriction rstr_j;
1549
1550 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j));
1551 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j));
1552 if (rstr_i == rstr_j && vec_i == vec_j) {
1553 field_rstr_in_buffer[j] = i;
1554 is_ordered[j] = true;
1555 input_field_order[curr_index] = j;
1556 curr_index++;
1557 }
1558 CeedCallBackend(CeedVectorDestroy(&vec_j));
1559 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
1560 }
1561 CeedCallBackend(CeedVectorDestroy(&vec_i));
1562 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
1563 }
1564 }
1565
1566 // -- Input restriction and basis
1567 code << "\n" << tab << "// -- Input field restrictions and basis actions\n";
1568 for (CeedInt i = 0; i < num_input_fields; i++) {
1569 const char *field_name;
1570 const CeedInt f = input_field_order[i];
1571
1572 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name));
1573 code << tab << "// ---- Input field " << f << ": " << field_name << "\n";
1574
1575 // ---- Restriction
1576 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f],
1577 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices));
1578
1579 // ---- Basis action
1580 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true,
1581 is_all_tensor, is_at_points, use_3d_slices));
1582 }
1583
1584 // -- Q function
1585 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
1586 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
1587 Q_1d, is_all_tensor, is_at_points, use_3d_slices, false));
1588
1589 // -- Output basis and restriction
1590 code << "\n" << tab << "// -- Output field basis action and restrictions\n";
1591 for (CeedInt i = 0; i < num_output_fields; i++) {
1592 const char *field_name;
1593
1594 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
1595 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
1596
1597 // ---- Basis action
1598 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d, false,
1599 is_all_tensor, is_at_points, use_3d_slices));
1600
1601 // ---- Restriction
1602 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, i, NULL, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d,
1603 false, is_all_tensor, is_at_points, use_3d_slices));
1604 }
1605
1606 // Close loop and function
1607 tab.pop();
1608 code << tab << "}\n";
1609 tab.pop();
1610 code << tab << "}\n";
1611 code << tab << "// -----------------------------------------------------------------------------\n\n";
1612
1613 // Compile
1614 {
1615 bool is_compile_good = false;
1616 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d);
1617
1618 data->thread_1d = T_1d;
1619 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good, &data->module, 1, "OP_T_1D", T_1d));
1620 if (is_compile_good) {
1621 *is_good_build = true;
1622 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module, operator_name.c_str(), &data->op));
1623 } else {
1624 *is_good_build = false;
1625 data->use_fallback = true;
1626 }
1627 }
1628 CeedCallBackend(CeedOperatorSetSetupDone(op));
1629 CeedCallBackend(CeedDestroy(&ceed));
1630 CeedCallBackend(CeedQFunctionDestroy(&qf));
1631 return CEED_ERROR_SUCCESS;
1632 }
1633
1634 //------------------------------------------------------------------------------
1635 // Build AtPoints assembly operator kernel
1636 //------------------------------------------------------------------------------
CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op,bool is_full,bool * is_good_build)1637 static int CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(CeedOperator op, bool is_full, bool *is_good_build) {
1638 bool is_all_tensor = true, is_at_points = false, use_3d_slices = false;
1639 Ceed ceed;
1640 CeedInt Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0, coords_comp_stride = 0;
1641 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1642 CeedQFunction_Cuda_gen *qf_data;
1643 CeedQFunction qf;
1644 CeedOperatorField *op_input_fields, *op_output_fields;
1645 CeedOperator_Cuda_gen *data;
1646 std::ostringstream code;
1647 Tab tab;
1648
1649 // Check compatibility
1650 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1651 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points));
1652 CeedCheck(is_at_points, ceed, CEED_ERROR_BACKEND, "Only AtPoints operator assembly supported");
1653
1654 // Retrieve operator data
1655 CeedCallBackend(CeedOperatorGetData(op, &data));
1656 Q = data->Q;
1657 Q_1d = data->Q_1d;
1658 max_dim = data->dim;
1659 {
1660 CeedElemRestriction rstr_points = NULL;
1661
1662 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
1663 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
1664 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr_points, &coords_comp_stride));
1665 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1666 }
1667 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1668 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data));
1669 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1670 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1671
1672 // Add atomicAdd function for old NVidia architectures
1673 {
1674 Ceed_Cuda *ceed_data;
1675 struct cudaDeviceProp prop;
1676
1677 CeedCallBackend(CeedGetData(ceed, &ceed_data));
1678 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id));
1679 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) {
1680 code << tab << "// AtomicAdd fallback source\n";
1681 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n";
1682 }
1683 }
1684
1685 // Load basis source files
1686 code << tab << "// Tensor basis source\n";
1687 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n";
1688 code << tab << "// AtPoints basis source\n";
1689 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h>\n\n";
1690 code << tab << "// CodeGen operator source\n";
1691 code << tab << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n";
1692
1693 // Get QFunction name
1694 std::string qfunction_name(qf_data->qfunction_name);
1695 std::string operator_name;
1696
1697 if (is_full) {
1698 operator_name = "CeedKernelCudaGenOperatorFullAssembly_" + qfunction_name;
1699 } else {
1700 operator_name = "CeedKernelCudaGenOperatorDiagonalAssembly_" + qfunction_name;
1701 }
1702
1703 // Define CEED_Q_VLA
1704 code << "\n" << tab << "#undef CEED_Q_VLA\n";
1705 code << tab << "#define CEED_Q_VLA 1\n\n";
1706
1707 // Add user QFunction source
1708 {
1709 const char *source_path;
1710
1711 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path));
1712 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file");
1713
1714 code << tab << "// User QFunction source\n";
1715 code << tab << "#include \"" << source_path << "\"\n\n";
1716 }
1717
1718 // Setup
1719 code << "\n" << tab << "// -----------------------------------------------------------------------------\n";
1720 code << tab << "// Operator Assembly Kernel\n";
1721 code << tab << "// \n";
1722 code << tab << "// d_[in,out]_i: CeedVector device array\n";
1723 code << tab << "// r_[in,out]_e_i: Element vector register\n";
1724 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n";
1725 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n";
1726 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n";
1727 code << tab << "// \n";
1728 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n";
1729 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n";
1730 code << tab << "// -----------------------------------------------------------------------------\n";
1731 code << tab << "extern \"C\" __global__ void " << operator_name
1732 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda "
1733 "points, CeedScalar *__restrict__ values_array) {\n";
1734 tab.push();
1735
1736 // Scratch buffers
1737 for (CeedInt i = 0; i < num_input_fields; i++) {
1738 CeedEvalMode eval_mode;
1739
1740 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1741 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT
1742 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n";
1743 }
1744 }
1745 for (CeedInt i = 0; i < num_output_fields; i++) {
1746 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n";
1747 }
1748
1749 code << tab << "const CeedInt max_dim = " << max_dim << ";\n";
1750 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n";
1751 code << tab << "const CeedInt max_num_points = " << max_num_points << ";\n";
1752 code << tab << "const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n";
1753
1754 // Shared data
1755 code << tab << "extern __shared__ CeedScalar slice[];\n";
1756 code << tab << "SharedData_Cuda data;\n";
1757 code << tab << "data.t_id_x = threadIdx.x;\n";
1758 code << tab << "data.t_id_y = threadIdx.y;\n";
1759 code << tab << "data.t_id_z = threadIdx.z;\n";
1760 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n";
1761 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n";
1762
1763 // -- Determine input mat reuse
1764 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX];
1765
1766 for (CeedInt i = 0; i < num_input_fields; i++) {
1767 input_matrix_reuse[i].index = -1;
1768 }
1769 for (CeedInt i = 0; i < num_input_fields; i++) {
1770 CeedEvalMode eval_mode_i;
1771 CeedBasis basis_i;
1772
1773 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i));
1774 if (eval_mode_i == CEED_EVAL_WEIGHT) continue;
1775 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i));
1776 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) {
1777 CeedEvalMode eval_mode_j;
1778 CeedBasis basis_j;
1779
1780 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
1781 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1782 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
1783 if (basis_i == basis_j) {
1784 input_matrix_reuse[i].index = j;
1785 input_matrix_reuse[i].is_input = true;
1786 input_matrix_reuse[i].eval_mode = eval_mode_j;
1787 }
1788 CeedCallBackend(CeedBasisDestroy(&basis_j));
1789 }
1790 CeedCallBackend(CeedBasisDestroy(&basis_i));
1791 }
1792
1793 // -- Determine output mat reuse
1794 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX];
1795
1796 for (CeedInt i = 0; i < num_output_fields; i++) {
1797 output_matrix_reuse[i].index = -1;
1798 }
1799 for (CeedInt i = 0; i < num_output_fields; i++) {
1800 CeedEvalMode eval_mode_i;
1801 CeedBasis basis_i;
1802
1803 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i));
1804 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i));
1805 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) {
1806 CeedEvalMode eval_mode_j;
1807 CeedBasis basis_j;
1808
1809 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
1810 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1811 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
1812 if (basis_i == basis_j) {
1813 output_matrix_reuse[i].index = j;
1814 output_matrix_reuse[i].is_input = true;
1815 output_matrix_reuse[i].eval_mode = eval_mode_j;
1816 }
1817 CeedCallBackend(CeedBasisDestroy(&basis_j));
1818 }
1819 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) {
1820 CeedEvalMode eval_mode_j;
1821 CeedBasis basis_j;
1822
1823 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j));
1824 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
1825 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j));
1826 if (basis_i == basis_j) {
1827 output_matrix_reuse[i].index = j;
1828 output_matrix_reuse[i].is_input = false;
1829 output_matrix_reuse[i].eval_mode = eval_mode_j;
1830 }
1831 CeedCallBackend(CeedBasisDestroy(&basis_j));
1832 }
1833 CeedCallBackend(CeedBasisDestroy(&basis_i));
1834 }
1835
1836 // Initialize constants, and matrices B and G
1837 code << "\n" << tab << "// Input field constants and basis data\n";
1838 for (CeedInt i = 0; i < num_input_fields; i++) {
1839 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i],
1840 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, false));
1841 }
1842 code << "\n" << tab << "// Output field constants and basis data\n";
1843 for (CeedInt i = 0; i < num_output_fields; i++) {
1844 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i],
1845 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, false));
1846 }
1847
1848 // Loop over all elements
1849 code << "\n" << tab << "// Element loop\n";
1850 code << tab << "__syncthreads();\n";
1851 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n";
1852 tab.push();
1853
1854 // -- Compute minimum buffer space needed
1855 CeedInt max_rstr_buffer_size = 1;
1856
1857 for (CeedInt i = 0; i < num_input_fields; i++) {
1858 CeedEvalMode eval_mode;
1859
1860 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
1861 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) {
1862 CeedInt num_comp;
1863 CeedElemRestriction elem_rstr;
1864
1865 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1866 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1867 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
1868 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1869 }
1870 }
1871 for (CeedInt i = 0; i < num_output_fields; i++) {
1872 CeedEvalMode eval_mode;
1873
1874 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
1875 if (eval_mode != CEED_EVAL_NONE) {
1876 CeedInt num_comp;
1877 CeedElemRestriction elem_rstr;
1878
1879 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
1880 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1881 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
1882 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1883 }
1884 }
1885 code << tab << "// Scratch restriction buffer space\n";
1886 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n";
1887
1888 // -- Determine best input field processing order
1889 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX];
1890
1891 for (CeedInt i = 0; i < num_input_fields; i++) {
1892 field_rstr_in_buffer[i] = -1;
1893 input_field_order[i] = -1;
1894 }
1895 {
1896 bool is_ordered[CEED_FIELD_MAX];
1897 CeedInt curr_index = 0;
1898
1899 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false;
1900 for (CeedInt i = 0; i < num_input_fields; i++) {
1901 CeedVector vec_i;
1902 CeedElemRestriction rstr_i;
1903
1904 if (is_ordered[i]) continue;
1905 field_rstr_in_buffer[i] = i;
1906 is_ordered[i] = true;
1907 input_field_order[curr_index] = i;
1908 curr_index++;
1909 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i));
1910 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT
1911 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i));
1912 for (CeedInt j = i + 1; j < num_input_fields; j++) {
1913 CeedVector vec_j;
1914 CeedElemRestriction rstr_j;
1915
1916 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j));
1917 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j));
1918 if (rstr_i == rstr_j && vec_i == vec_j) {
1919 field_rstr_in_buffer[j] = i;
1920 is_ordered[j] = true;
1921 input_field_order[curr_index] = j;
1922 curr_index++;
1923 }
1924 CeedCallBackend(CeedVectorDestroy(&vec_j));
1925 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
1926 }
1927 CeedCallBackend(CeedVectorDestroy(&vec_i));
1928 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
1929 }
1930 }
1931
1932 // -- Input restriction and basis
1933 code << "\n" << tab << "// -- Input field restrictions and basis actions\n";
1934 CeedInt active_field_index = -1;
1935
1936 for (CeedInt i = 0; i < num_input_fields; i++) {
1937 bool is_active = false;
1938 const char *field_name;
1939 const CeedInt f = input_field_order[i];
1940
1941 {
1942 CeedVector vec;
1943
1944 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec));
1945 is_active = vec == CEED_VECTOR_ACTIVE;
1946 CeedCallBackend(CeedVectorDestroy(&vec));
1947 }
1948
1949 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name));
1950 code << tab << "// ---- Input field " << f << ": " << field_name << "\n";
1951
1952 if (is_active) {
1953 std::string var_suffix = "_in_" + std::to_string(f);
1954
1955 code << tab << "// Active field - no restriction or basis action here\n";
1956 if (active_field_index == -1) {
1957 active_field_index = f;
1958 code << tab << "CeedScalar r_e" << var_suffix << "[num_comp" << var_suffix << "*" << (max_dim >= 3 ? "P_1d" + var_suffix : "1")
1959 << "] = {0.0};\n";
1960 } else {
1961 code << tab << "CeedScalar *r_e" << var_suffix << " = r_e_in_" << active_field_index << ";\n";
1962 }
1963 } else {
1964 // ---- Restriction
1965 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f],
1966 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices));
1967
1968 // ---- Basis action
1969 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true,
1970 is_all_tensor, is_at_points, use_3d_slices));
1971 }
1972 }
1973
1974 // -- Loop over active field
1975 std::string active_var_suffix = "_in_" + std::to_string(active_field_index);
1976
1977 code << "\n" << tab << "// Loop over nodes in active field\n";
1978 code << tab << "for (CeedInt n = 0; n < num_comp" << active_var_suffix << "*P_1d" << active_var_suffix
1979 << (max_dim > 1 ? "*P_1d" + active_var_suffix : "") << (max_dim > 2 ? "*P_1d" + active_var_suffix : "") << "; n++) {\n";
1980 tab.push();
1981
1982 // -- Set current active node and component to 1
1983 code << tab << "// Set current active node and component to 1.0\n";
1984 code << tab << "SetEVecStandard" << max_dim << "d_Single<num_comp" << active_var_suffix << ", P_1d" << active_var_suffix << ">(data, n, 1.0, r_e"
1985 << active_var_suffix << ");\n\n";
1986
1987 for (CeedInt i = 0; i < num_input_fields; i++) {
1988 bool is_active = false;
1989 const char *field_name;
1990 const CeedInt f = input_field_order[i];
1991
1992 {
1993 CeedVector vec;
1994
1995 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec));
1996 is_active = vec == CEED_VECTOR_ACTIVE;
1997 CeedCallBackend(CeedVectorDestroy(&vec));
1998 }
1999 if (!is_active) continue;
2000
2001 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name));
2002 code << tab << "// ---- Input field " << f << ": " << field_name << "\n";
2003
2004 // ---- Basis action
2005 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true,
2006 is_all_tensor, is_at_points, use_3d_slices));
2007 }
2008
2009 // -- Q function
2010 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
2011 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
2012 Q_1d, is_all_tensor, is_at_points, use_3d_slices, true));
2013
2014 // -- Output basis and restriction
2015 code << "\n" << tab << "// -- Output field basis action and restrictions\n";
2016 for (CeedInt i = 0; i < num_output_fields; i++) {
2017 bool is_active = false;
2018 const char *field_name;
2019
2020 {
2021 CeedVector vec;
2022
2023 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
2024 is_active = vec == CEED_VECTOR_ACTIVE;
2025 CeedCallBackend(CeedVectorDestroy(&vec));
2026 }
2027 if (!is_active) continue;
2028
2029 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
2030 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
2031
2032 // ---- Basis action
2033 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], max_dim, Q_1d, false,
2034 is_all_tensor, is_at_points, use_3d_slices));
2035
2036 // ---- Restriction
2037 if (is_full) {
2038 std::string var_suffix = "_out_" + std::to_string(i);
2039 CeedInt comp_stride;
2040 CeedSize l_size;
2041 CeedElemRestriction elem_rstr;
2042
2043 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
2044 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size));
2045 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n";
2046 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
2047 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
2048 code << tab << "WriteLVecStandard" << max_dim << "d_Assembly<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", P_1d" + var_suffix
2049 << ">(data, l_size" << var_suffix << ", elem, n, r_e" << var_suffix << ", values_array);\n";
2050 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
2051 } else {
2052 std::string var_suffix = "_out_" + std::to_string(i);
2053 CeedInt comp_stride;
2054 CeedSize l_size;
2055 CeedElemRestriction elem_rstr;
2056
2057 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
2058 CeedCallBackend(CeedElemRestrictionGetLVectorSize(elem_rstr, &l_size));
2059 code << tab << "const CeedInt l_size" << var_suffix << " = " << l_size << ";\n";
2060 CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride));
2061 code << tab << "const CeedInt comp_stride" << var_suffix << " = " << comp_stride << ";\n";
2062 code << tab << "WriteLVecStandard" << max_dim << "d_Single<num_comp" << var_suffix << ", comp_stride" << var_suffix << ", P_1d" + var_suffix
2063 << ">(data, l_size" << var_suffix << ", elem, n, indices.outputs[" << i << "], r_e" << var_suffix << ", values_array);\n";
2064 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
2065 }
2066 }
2067
2068 // -- Reset current active node and component
2069 code << "\n" << tab << "// Reset current active node and component to 0.0\n";
2070 code << tab << "SetEVecStandard" << max_dim << "d_Single<num_comp" << active_var_suffix << ", P_1d" << active_var_suffix << ">(data, n, 0.0, r_e"
2071 << active_var_suffix << ");\n";
2072
2073 // -- End of loop over active field
2074 tab.pop();
2075 code << tab << "}\n";
2076
2077 // Close loop and function
2078 tab.pop();
2079 code << tab << "}\n";
2080 tab.pop();
2081 code << tab << "}\n";
2082 code << tab << "// -----------------------------------------------------------------------------\n\n";
2083
2084 // Compile
2085 {
2086 bool is_compile_good = false;
2087 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d);
2088
2089 data->thread_1d = T_1d;
2090 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good,
2091 is_full ? &data->module_assemble_full : &data->module_assemble_diagonal, 1, "OP_T_1D", T_1d));
2092 if (is_compile_good) {
2093 *is_good_build = true;
2094 CeedCallBackend(CeedGetKernel_Cuda(ceed, is_full ? data->module_assemble_full : data->module_assemble_diagonal, operator_name.c_str(),
2095 is_full ? &data->assemble_full : &data->assemble_diagonal));
2096 } else {
2097 *is_good_build = false;
2098 data->use_assembly_fallback = true;
2099 }
2100 }
2101 CeedCallBackend(CeedDestroy(&ceed));
2102 CeedCallBackend(CeedQFunctionDestroy(&qf));
2103 return CEED_ERROR_SUCCESS;
2104 }
2105
CeedOperatorBuildKernelDiagonalAssemblyAtPoints_Cuda_gen(CeedOperator op,bool * is_good_build)2106 extern "C" int CeedOperatorBuildKernelDiagonalAssemblyAtPoints_Cuda_gen(CeedOperator op, bool *is_good_build) {
2107 return CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(op, false, is_good_build);
2108 }
2109
CeedOperatorBuildKernelFullAssemblyAtPoints_Cuda_gen(CeedOperator op,bool * is_good_build)2110 extern "C" int CeedOperatorBuildKernelFullAssemblyAtPoints_Cuda_gen(CeedOperator op, bool *is_good_build) {
2111 return CeedOperatorBuildKernelAssemblyAtPoints_Cuda_gen(op, true, is_good_build);
2112 }
2113
2114 //------------------------------------------------------------------------------
2115 // Build QFunction assembly operator kernel
2116 //------------------------------------------------------------------------------
CeedOperatorBuildKernelLinearAssembleQFunction_Cuda_gen(CeedOperator op,bool * is_good_build)2117 extern "C" int CeedOperatorBuildKernelLinearAssembleQFunction_Cuda_gen(CeedOperator op, bool *is_good_build) {
2118 bool is_all_tensor = true, is_all_nontensor = true, is_at_points = false, use_3d_slices = false;
2119 Ceed ceed;
2120 CeedInt Q, Q_1d, num_input_fields, num_output_fields, max_dim = 1, max_num_points = 0;
2121 CeedQFunctionField *qf_input_fields, *qf_output_fields;
2122 CeedQFunction_Cuda_gen *qf_data;
2123 CeedQFunction qf;
2124 CeedOperatorField *op_input_fields, *op_output_fields;
2125 CeedOperator_Cuda_gen *data;
2126 std::ostringstream code;
2127 Tab tab;
2128
2129 // Check compatibility
2130 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
2131 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points));
2132 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "AtPoints QFunction assembly is not supported");
2133
2134 // Check field compatibility
2135 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
2136 {
2137 bool has_shared_bases = true;
2138
2139 for (CeedInt i = 0; i < num_input_fields; i++) {
2140 CeedBasis basis;
2141
2142 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
2143 if (basis != CEED_BASIS_NONE) {
2144 bool is_tensor = true;
2145 const char *resource;
2146 char *resource_root;
2147 Ceed basis_ceed;
2148
2149 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
2150 is_all_tensor = is_all_tensor && is_tensor;
2151 is_all_nontensor = is_all_nontensor && !is_tensor;
2152 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed));
2153 CeedCallBackend(CeedGetResource(basis_ceed, &resource));
2154 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root));
2155 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared");
2156 CeedCallBackend(CeedFree(&resource_root));
2157 CeedCallBackend(CeedDestroy(&basis_ceed));
2158 }
2159 CeedCallBackend(CeedBasisDestroy(&basis));
2160 }
2161
2162 for (CeedInt i = 0; i < num_output_fields; i++) {
2163 CeedBasis basis;
2164
2165 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
2166 if (basis != CEED_BASIS_NONE) {
2167 bool is_tensor = true;
2168 const char *resource;
2169 char *resource_root;
2170 Ceed basis_ceed;
2171
2172 CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor));
2173 is_all_tensor = is_all_tensor && is_tensor;
2174 is_all_nontensor = is_all_nontensor && !is_tensor;
2175
2176 CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed));
2177 CeedCallBackend(CeedGetResource(basis_ceed, &resource));
2178 CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root));
2179 has_shared_bases = has_shared_bases && !strcmp(resource_root, "/gpu/cuda/shared");
2180 CeedCallBackend(CeedFree(&resource_root));
2181 CeedCallBackend(CeedDestroy(&basis_ceed));
2182 }
2183 CeedCallBackend(CeedBasisDestroy(&basis));
2184 }
2185 }
2186
2187 // Retrieve operator data
2188 CeedCallBackend(CeedOperatorGetData(op, &data));
2189 Q = data->Q;
2190 Q_1d = data->Q_1d;
2191 max_dim = data->dim;
2192 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
2193 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data));
2194 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
2195
2196 // Add atomicAdd function for old NVidia architectures
2197 {
2198 Ceed_Cuda *ceed_data;
2199 struct cudaDeviceProp prop;
2200
2201 CeedCallBackend(CeedGetData(ceed, &ceed_data));
2202 CeedCallBackend(cudaGetDeviceProperties(&prop, ceed_data->device_id));
2203 if ((prop.major < 6) && (CEED_SCALAR_TYPE != CEED_SCALAR_FP32)) {
2204 code << tab << "// AtomicAdd fallback source\n";
2205 code << tab << "#include <ceed/jit-source/cuda/cuda-atomic-add-fallback.h>\n\n";
2206 }
2207 }
2208
2209 // Load basis source files
2210 if (!is_all_nontensor) {
2211 code << tab << "// Tensor basis source\n";
2212 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-templates.h>\n\n";
2213 }
2214 if (!is_all_tensor) {
2215 code << tab << "// Non-tensor basis source\n";
2216 code << tab << "#include <ceed/jit-source/cuda/cuda-shared-basis-nontensor-templates.h>\n\n";
2217 }
2218 if (!is_all_tensor && !is_all_nontensor) {
2219 code << "// Tensor basis source\n";
2220 code << "#include <ceed/jit-source/cuda/cuda-shared-basis-tensor-flattened-templates.h>\n\n";
2221 }
2222 code << "// CodeGen operator source\n";
2223 code << "#include <ceed/jit-source/cuda/cuda-gen-templates.h>\n\n";
2224
2225 // Get QFunction name
2226 std::string qfunction_name(qf_data->qfunction_name);
2227 std::string operator_name;
2228
2229 operator_name = "CeedKernelCudaGenQFunctionAssembly_" + qfunction_name;
2230
2231 // Define CEED_Q_VLA
2232 code << "\n" << tab << "#undef CEED_Q_VLA\n";
2233 if (max_dim != 3 || is_at_points || use_3d_slices || !is_all_tensor) {
2234 code << tab << "#define CEED_Q_VLA 1\n\n";
2235 } else {
2236 code << tab << "#define CEED_Q_VLA " << Q_1d << "\n\n";
2237 }
2238
2239 // Add user QFunction source
2240 {
2241 const char *source_path;
2242
2243 CeedCallBackend(CeedQFunctionGetSourcePath(qf, &source_path));
2244 CeedCheck(source_path, ceed, CEED_ERROR_UNSUPPORTED, "/gpu/cuda/gen backend requires QFunction source code file");
2245
2246 code << tab << "// User QFunction source\n";
2247 code << tab << "#include \"" << source_path << "\"\n\n";
2248 }
2249
2250 // Setup
2251 code << "\n" << tab << "// -----------------------------------------------------------------------------\n";
2252 code << tab << "// Operator Assembly Kernel\n";
2253 code << tab << "// \n";
2254 code << tab << "// d_[in,out]_i: CeedVector device array\n";
2255 code << tab << "// r_[in,out]_e_i: Element vector register\n";
2256 code << tab << "// r_[in,out]_q_i: Quadrature space vector register\n";
2257 code << tab << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n";
2258 code << tab << "// r_[in,out]_s_i: Quadrature space slice vector register\n";
2259 code << tab << "// \n";
2260 code << tab << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n";
2261 code << tab << "// s_G_[in,out]_i: Gradient matrix, shared memory\n";
2262 code << tab << "// -----------------------------------------------------------------------------\n";
2263 code << tab << "extern \"C\" __global__ void " << operator_name
2264 << "(CeedInt num_elem, void* ctx, FieldsInt_Cuda indices, Fields_Cuda fields, Fields_Cuda B, Fields_Cuda G, CeedScalar *W, Points_Cuda "
2265 "points, CeedScalar *__restrict__ values_array) {\n";
2266 tab.push();
2267
2268 // Scratch buffers
2269 for (CeedInt i = 0; i < num_input_fields; i++) {
2270 CeedEvalMode eval_mode;
2271
2272 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
2273 if (eval_mode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT
2274 code << tab << "const CeedScalar *__restrict__ d_in_" << i << " = fields.inputs[" << i << "];\n";
2275 }
2276 }
2277 for (CeedInt i = 0; i < num_output_fields; i++) {
2278 bool is_active = false;
2279
2280 {
2281 CeedVector vec;
2282
2283 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
2284 is_active = vec == CEED_VECTOR_ACTIVE;
2285 CeedCallBackend(CeedVectorDestroy(&vec));
2286 }
2287 if (is_active) {
2288 code << tab << "CeedScalar *__restrict__ d_out_" << i << " = fields.outputs[" << i << "];\n";
2289 }
2290 }
2291
2292 code << tab << "const CeedInt max_dim = " << max_dim << ";\n";
2293 if (!is_all_tensor) {
2294 code << tab << "const CeedInt Q = " << Q << ";\n";
2295 }
2296 if (!is_all_nontensor) {
2297 code << tab << "const CeedInt Q_1d = " << Q_1d << ";\n";
2298 }
2299
2300 // Shared data
2301 code << tab << "extern __shared__ CeedScalar slice[];\n";
2302 code << tab << "SharedData_Cuda data;\n";
2303 code << tab << "data.t_id_x = threadIdx.x;\n";
2304 code << tab << "data.t_id_y = threadIdx.y;\n";
2305 code << tab << "data.t_id_z = threadIdx.z;\n";
2306 code << tab << "data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n";
2307 code << tab << "data.slice = slice + data.t_id_z*OP_T_1D" << ((!is_all_tensor || max_dim == 1) ? "" : "*OP_T_1D") << ";\n";
2308
2309 // -- Determine input mat reuse
2310 FieldReuse_Cuda input_matrix_reuse[CEED_FIELD_MAX];
2311
2312 for (CeedInt i = 0; i < num_input_fields; i++) {
2313 input_matrix_reuse[i].index = -1;
2314 }
2315 for (CeedInt i = 0; i < num_input_fields; i++) {
2316 bool is_tensor = true;
2317 CeedEvalMode eval_mode_i;
2318 CeedBasis basis_i;
2319
2320 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode_i));
2321 if (eval_mode_i == CEED_EVAL_WEIGHT) continue;
2322 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis_i));
2323 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor));
2324 for (CeedInt j = 0; (input_matrix_reuse[i].index == -1) && (j < i); j++) {
2325 CeedEvalMode eval_mode_j;
2326 CeedBasis basis_j;
2327
2328 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
2329 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
2330 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
2331 if (basis_i == basis_j) {
2332 if (is_tensor) {
2333 input_matrix_reuse[i].index = j;
2334 input_matrix_reuse[i].is_input = true;
2335 input_matrix_reuse[i].eval_mode = eval_mode_j;
2336 } else {
2337 // For non-tensor can only re-use with the same eval mode
2338 if (eval_mode_i == eval_mode_j) {
2339 input_matrix_reuse[i].index = j;
2340 input_matrix_reuse[i].is_input = true;
2341 input_matrix_reuse[i].eval_mode = eval_mode_j;
2342 }
2343 }
2344 }
2345 CeedCallBackend(CeedBasisDestroy(&basis_j));
2346 }
2347 CeedCallBackend(CeedBasisDestroy(&basis_i));
2348 }
2349
2350 // -- Determine output mat reuse
2351 FieldReuse_Cuda output_matrix_reuse[CEED_FIELD_MAX];
2352
2353 for (CeedInt i = 0; i < num_output_fields; i++) {
2354 output_matrix_reuse[i].index = -1;
2355 }
2356 for (CeedInt i = 0; i < num_output_fields; i++) {
2357 bool is_tensor = true;
2358 CeedEvalMode eval_mode_i;
2359 CeedBasis basis_i;
2360
2361 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode_i));
2362 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis_i));
2363 CeedCallBackend(CeedBasisIsTensor(basis_i, &is_tensor));
2364 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < num_input_fields); j++) {
2365 CeedEvalMode eval_mode_j;
2366 CeedBasis basis_j;
2367
2368 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[j], &eval_mode_j));
2369 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
2370 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[j], &basis_j));
2371 if (basis_i == basis_j) {
2372 if (is_tensor) {
2373 output_matrix_reuse[i].index = j;
2374 output_matrix_reuse[i].is_input = true;
2375 output_matrix_reuse[i].eval_mode = eval_mode_j;
2376 } else {
2377 // For non-tensor can only re-use with the same eval mode
2378 if (eval_mode_i == eval_mode_j) {
2379 output_matrix_reuse[i].index = j;
2380 output_matrix_reuse[i].is_input = true;
2381 output_matrix_reuse[i].eval_mode = eval_mode_j;
2382 }
2383 }
2384 }
2385 CeedCallBackend(CeedBasisDestroy(&basis_j));
2386 }
2387 for (CeedInt j = 0; (output_matrix_reuse[i].index == -1) && (j < i); j++) {
2388 CeedEvalMode eval_mode_j;
2389 CeedBasis basis_j;
2390
2391 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[j], &eval_mode_j));
2392 if (eval_mode_j == CEED_EVAL_WEIGHT) continue;
2393 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[j], &basis_j));
2394 if (basis_i == basis_j) {
2395 if (is_tensor) {
2396 output_matrix_reuse[i].index = j;
2397 output_matrix_reuse[i].is_input = false;
2398 output_matrix_reuse[i].eval_mode = eval_mode_j;
2399 } else {
2400 // For non-tensor can only re-use with the same eval mode
2401 if (eval_mode_i == eval_mode_j) {
2402 output_matrix_reuse[i].index = j;
2403 output_matrix_reuse[i].is_input = false;
2404 output_matrix_reuse[i].eval_mode = eval_mode_j;
2405 }
2406 }
2407 }
2408 CeedCallBackend(CeedBasisDestroy(&basis_j));
2409 }
2410 CeedCallBackend(CeedBasisDestroy(&basis_i));
2411 }
2412
2413 // Initialize constants, and matrices B and G
2414 code << "\n" << tab << "// Input field constants and basis data\n";
2415 for (CeedInt i = 0; i < num_input_fields; i++) {
2416 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_input_fields[i], qf_input_fields[i], input_matrix_reuse[i],
2417 max_dim, Q, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices, true));
2418 }
2419 code << "\n" << tab << "// Output field constants and basis data\n";
2420 for (CeedInt i = 0; i < num_output_fields; i++) {
2421 CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, tab, i, op_output_fields[i], qf_output_fields[i], output_matrix_reuse[i],
2422 max_dim, Q, Q_1d, false, is_all_tensor, is_at_points, use_3d_slices, true));
2423 }
2424
2425 // Loop over all elements
2426 code << "\n" << tab << "// Element loop\n";
2427 code << tab << "__syncthreads();\n";
2428 code << tab << "for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < num_elem; elem += gridDim.x*blockDim.z) {\n";
2429 tab.push();
2430
2431 // -- Compute minimum buffer space needed
2432 CeedInt max_rstr_buffer_size = 1;
2433
2434 for (CeedInt i = 0; i < num_input_fields; i++) {
2435 CeedEvalMode eval_mode;
2436
2437 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
2438 if (eval_mode != CEED_EVAL_NONE && eval_mode != CEED_EVAL_WEIGHT) {
2439 CeedInt num_comp;
2440 CeedElemRestriction elem_rstr;
2441
2442 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
2443 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
2444 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
2445 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
2446 }
2447 }
2448 for (CeedInt i = 0; i < num_output_fields; i++) {
2449 CeedEvalMode eval_mode;
2450
2451 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
2452 if (eval_mode != CEED_EVAL_NONE) {
2453 CeedInt num_comp;
2454 CeedElemRestriction elem_rstr;
2455
2456 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
2457 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
2458 max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_all_tensor && (max_dim >= 3) ? Q_1d : 1));
2459 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
2460 }
2461 }
2462 code << tab << "// Scratch restriction buffer space\n";
2463 code << tab << "CeedScalar r_e_scratch[" << max_rstr_buffer_size << "];\n";
2464
2465 // -- Determine best input field processing order
2466 CeedInt field_rstr_in_buffer[CEED_FIELD_MAX], input_field_order[CEED_FIELD_MAX];
2467
2468 for (CeedInt i = 0; i < num_input_fields; i++) {
2469 field_rstr_in_buffer[i] = -1;
2470 input_field_order[i] = -1;
2471 }
2472 {
2473 bool is_ordered[CEED_FIELD_MAX];
2474 CeedInt curr_index = 0;
2475
2476 for (CeedInt i = 0; i < num_input_fields; i++) is_ordered[i] = false;
2477 for (CeedInt i = 0; i < num_input_fields; i++) {
2478 CeedVector vec_i;
2479 CeedElemRestriction rstr_i;
2480
2481 if (is_ordered[i]) continue;
2482 field_rstr_in_buffer[i] = i;
2483 is_ordered[i] = true;
2484 input_field_order[curr_index] = i;
2485 curr_index++;
2486 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec_i));
2487 if (vec_i == CEED_VECTOR_NONE) continue; // CEED_EVAL_WEIGHT
2488 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &rstr_i));
2489 for (CeedInt j = i + 1; j < num_input_fields; j++) {
2490 CeedVector vec_j;
2491 CeedElemRestriction rstr_j;
2492
2493 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec_j));
2494 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[j], &rstr_j));
2495 if (rstr_i == rstr_j && vec_i == vec_j) {
2496 field_rstr_in_buffer[j] = i;
2497 is_ordered[j] = true;
2498 input_field_order[curr_index] = j;
2499 curr_index++;
2500 }
2501 CeedCallBackend(CeedVectorDestroy(&vec_j));
2502 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
2503 }
2504 CeedCallBackend(CeedVectorDestroy(&vec_i));
2505 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
2506 }
2507 }
2508
2509 // -- Input restriction and basis
2510 code << "\n" << tab << "// -- Input field restrictions and basis actions\n";
2511 CeedInt num_active_in = 0, num_active_out = 0, qf_assembly_size_out = 0;
2512 CeedInt active_fields_in[CEED_FIELD_MAX], active_fields_out[CEED_FIELD_MAX];
2513
2514 for (CeedInt i = 0; i < num_input_fields; i++) {
2515 bool is_active = false;
2516 const char *field_name;
2517 const CeedInt f = input_field_order[i];
2518
2519 {
2520 CeedVector vec;
2521
2522 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[f], &vec));
2523 is_active = vec == CEED_VECTOR_ACTIVE;
2524 CeedCallBackend(CeedVectorDestroy(&vec));
2525 }
2526
2527 CeedCallBackend(CeedOperatorFieldGetName(op_input_fields[f], &field_name));
2528 code << tab << "// ---- Input field " << f << ": " << field_name << "\n";
2529
2530 if (is_active) {
2531 CeedEvalMode eval_mode;
2532 CeedInt field_size;
2533
2534 active_fields_in[num_active_in] = f;
2535 num_active_in++;
2536 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[f], &field_size));
2537 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[f], &eval_mode));
2538 if (eval_mode == CEED_EVAL_GRAD) {
2539 code << tab << "CeedScalar r_q_in_" << f << "[num_comp_in_" << f << "*" << "dim_in_" << f << "*"
2540 << (is_all_tensor && (max_dim >= 3) ? "Q_1d" : "1") << "] = {0.};\n";
2541 } else {
2542 code << tab << "CeedScalar r_q_in_" << f << "[num_comp_in_" << f << "*" << (is_all_tensor && (max_dim >= 3) ? "Q_1d" : "1") << "] = {0.};\n";
2543 }
2544 code << tab << "const CeedInt field_size_in_" << f << " = " << field_size << ";\n";
2545 } else {
2546 // ---- Restriction
2547 CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, tab, f, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f],
2548 max_dim, Q_1d, true, is_all_tensor, is_at_points, use_3d_slices));
2549
2550 // ---- Basis action
2551 CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, tab, f, op_input_fields[f], qf_input_fields[f], max_dim, Q_1d, true,
2552 is_all_tensor, is_at_points, use_3d_slices));
2553 }
2554 }
2555 code << tab << "const CeedInt field_sizes_in[" << num_active_in << "] = {";
2556 for (CeedInt i = 0; i < num_active_in; i++) {
2557 code << "field_size_in_" << active_fields_in[i] << (i < num_active_in - 1 ? ", " : "");
2558 }
2559 code << "};\n";
2560 code << tab << "CeedScalar * r_q_in[" << num_active_in << "] = {";
2561 for (CeedInt i = 0; i < num_active_in; i++) {
2562 code << "r_q_in_" << active_fields_in[i] << (i < num_active_in - 1 ? ", " : "");
2563 }
2564 code << "};\n";
2565
2566 for (CeedInt i = 0; i < num_output_fields; i++) {
2567 bool is_active = false;
2568
2569 {
2570 CeedVector vec;
2571
2572 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
2573 is_active = vec == CEED_VECTOR_ACTIVE;
2574 CeedCallBackend(CeedVectorDestroy(&vec));
2575 }
2576 if (is_active) {
2577 const char *field_name;
2578 CeedInt field_size;
2579
2580 active_fields_out[num_active_out] = i;
2581 num_active_out++;
2582 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size));
2583 qf_assembly_size_out += field_size;
2584 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
2585 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
2586 code << tab << "const CeedInt field_size_out_" << i << " = " << field_size << ";\n";
2587 }
2588 }
2589 code << tab << "const CeedInt field_sizes_out[" << num_active_out << "] = {";
2590 for (CeedInt i = 0; i < num_active_out; i++) {
2591 code << "field_size_out_" << active_fields_out[i] << (i < num_active_out - 1 ? ", " : "");
2592 }
2593 code << "};\n";
2594 code << tab << "const CeedInt total_size_out = " << qf_assembly_size_out << ";\n";
2595
2596 // -- Loop over active field
2597 code << "\n" << tab << "CeedInt input_offset = 0;\n";
2598 code << tab << "// Loop over active QFunction input fields\n";
2599 code << tab << "const CeedInt num_active_in = " << num_active_in << ";\n";
2600 code << tab << "for (CeedInt a = 0; a < num_active_in; a++) {\n";
2601 tab.push();
2602
2603 // -- Loop over size of active field
2604 code << "\n" << tab << "// Loop over current active input field size\n";
2605 code << tab << "const CeedInt field_size_in = field_sizes_in[a];\n";
2606 code << tab << "for (CeedInt s = 0; s < field_size_in; s++) {\n";
2607 tab.push();
2608
2609 // -- Set current active point and component to 1
2610 code << tab << "// Set current active point and component to 1.0\n";
2611 if (is_all_tensor && (max_dim >= 3)) {
2612 code << tab << "for (CeedInt i = 0; i < Q_1d; i++) r_q_in[a][i + s * Q_1d] = 1.0;\n";
2613 } else {
2614 code << tab << "r_q_in[a][s] = 1.0;\n";
2615 }
2616
2617 // -- Q function
2618 CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, tab, max_dim, max_num_points, num_input_fields, op_input_fields,
2619 qf_input_fields, num_output_fields, op_output_fields, qf_output_fields, qfunction_name,
2620 Q_1d, is_all_tensor, is_at_points, use_3d_slices, true));
2621
2622 // -- Output basis and restriction
2623 code << "\n" << tab << "// -- Output field basis action and restrictions\n";
2624 CeedScalar offset = 0;
2625
2626 for (CeedInt i = 0; i < num_output_fields; i++) {
2627 bool is_active = false;
2628 const char *field_name;
2629
2630 {
2631 CeedVector vec;
2632
2633 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
2634 is_active = vec == CEED_VECTOR_ACTIVE;
2635 CeedCallBackend(CeedVectorDestroy(&vec));
2636 }
2637 if (!is_active) continue;
2638
2639 CeedCallBackend(CeedOperatorFieldGetName(op_output_fields[i], &field_name));
2640 code << tab << "// ---- Output field " << i << ": " << field_name << "\n";
2641
2642 // ---- Restriction
2643 CeedInt field_size;
2644
2645 code << tab << "WriteLVecStandard" << (is_all_tensor ? max_dim : 1) << "d_QFAssembly<total_size_out, field_size_out_" << i << ", "
2646 << (is_all_tensor ? "Q_1d" : "Q") << ">(data, num_elem, elem, input_offset + s, " << offset << ", r_q_out_" << i << ", values_array);\n";
2647 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size));
2648 offset += field_size;
2649 }
2650
2651 // -- Reset current active node and component
2652 code << "\n" << tab << "// Reset current active node and component to 0.0\n";
2653 if (is_all_tensor && (max_dim >= 3)) {
2654 code << tab << "for (CeedInt i = 0; i < Q_1d; i++) r_q_in[a][i + s * Q_1d] = 0.0;\n";
2655 } else {
2656 code << tab << "r_q_in[a][s] = 0.0;\n";
2657 }
2658
2659 // -- End of loop over size of active field
2660 tab.pop();
2661 code << tab << "}\n";
2662 code << tab << "input_offset += field_size_in;\n";
2663
2664 // -- End of loop over active field
2665 tab.pop();
2666 code << tab << "}\n";
2667
2668 // Close loop and function
2669 tab.pop();
2670 code << tab << "}\n";
2671 tab.pop();
2672 code << tab << "}\n";
2673 code << tab << "// -----------------------------------------------------------------------------\n\n";
2674
2675 // Compile
2676 {
2677 bool is_compile_good = false;
2678 const CeedInt T_1d = CeedIntMax(is_all_tensor ? Q_1d : Q, data->max_P_1d);
2679
2680 data->thread_1d = T_1d;
2681 CeedCallBackend(CeedTryCompile_Cuda(ceed, code.str().c_str(), &is_compile_good, &data->module_assemble_qfunction, 1, "OP_T_1D", T_1d));
2682 if (is_compile_good) {
2683 *is_good_build = true;
2684 CeedCallBackend(CeedGetKernel_Cuda(ceed, data->module_assemble_qfunction, operator_name.c_str(), &data->assemble_qfunction));
2685 } else {
2686 *is_good_build = false;
2687 data->use_assembly_fallback = true;
2688 }
2689 }
2690 CeedCallBackend(CeedDestroy(&ceed));
2691 CeedCallBackend(CeedQFunctionDestroy(&qf));
2692 return CEED_ERROR_SUCCESS;
2693 }
2694
2695 //------------------------------------------------------------------------------
2696