xref: /libCEED/backends/ref/ceed-ref-basis.c (revision 22ab0487938d6416bda03d32bba2b2245fabcc02) 
1 // Copyright (c) 2017-2024, 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 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <math.h>
11 #include <stdbool.h>
12 #include <string.h>
13 
14 #include "ceed-ref.h"
15 
16 //------------------------------------------------------------------------------
17 // Basis Apply
18 //------------------------------------------------------------------------------
19 static int CeedBasisApplyCore_Ref(CeedBasis basis, bool apply_add, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U,
20                                   CeedVector V) {
21   Ceed               ceed;
22   bool               is_tensor_basis, add = apply_add || (t_mode == CEED_TRANSPOSE);
23   CeedInt            dim, num_comp, q_comp, num_nodes, num_qpts;
24   const CeedScalar  *u;
25   CeedScalar        *v;
26   CeedTensorContract contract;
27   CeedBasis_Ref     *impl;
28 
29   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
30   CeedCallBackend(CeedBasisGetData(basis, &impl));
31   CeedCallBackend(CeedBasisGetDimension(basis, &dim));
32   CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
33   CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
34   CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes));
35   CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
36   CeedCallBackend(CeedBasisGetTensorContract(basis, &contract));
37   if (U != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u));
38   else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
39   // Clear v if operating in transpose
40   if (apply_add) CeedCallBackend(CeedVectorGetArray(V, CEED_MEM_HOST, &v));
41   else CeedCallBackend(CeedVectorGetArrayWrite(V, CEED_MEM_HOST, &v));
42 
43   if (t_mode == CEED_TRANSPOSE && !apply_add) {
44     CeedSize len;
45 
46     CeedCallBackend(CeedVectorGetLength(V, &len));
47     for (CeedInt i = 0; i < len; i++) v[i] = 0.0;
48   }
49 
50   CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor_basis));
51   if (is_tensor_basis) {
52     // Tensor basis
53     CeedInt P_1d, Q_1d;
54 
55     CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d));
56     CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
57     switch (eval_mode) {
58       // Interpolate to/from quadrature points
59       case CEED_EVAL_INTERP: {
60         if (impl->has_collo_interp) {
61           memcpy(v, u, num_elem * num_comp * num_nodes * sizeof(u[0]));
62         } else {
63           CeedInt P = P_1d, Q = Q_1d;
64 
65           if (t_mode == CEED_TRANSPOSE) {
66             P = Q_1d;
67             Q = P_1d;
68           }
69           CeedInt           pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
70           CeedScalar        tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)];
71           const CeedScalar *interp_1d;
72 
73           CeedCallBackend(CeedBasisGetInterp1D(basis, &interp_1d));
74           for (CeedInt d = 0; d < dim; d++) {
75             CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, interp_1d, t_mode, add && (d == dim - 1), d == 0 ? u : tmp[d % 2],
76                                                     d == dim - 1 ? v : tmp[(d + 1) % 2]));
77             pre /= P;
78             post *= Q;
79           }
80         }
81       } break;
82       // Evaluate the gradient to/from quadrature points
83       case CEED_EVAL_GRAD: {
84         // In CEED_NOTRANSPOSE mode:
85         // u has shape [dim, num_comp, P^dim, num_elem], row-major layout
86         // v has shape [dim, num_comp, Q^dim, num_elem], row-major layout
87         // In CEED_TRANSPOSE mode, the sizes of u and v are switched.
88         CeedInt P = P_1d, Q = Q_1d;
89 
90         if (t_mode == CEED_TRANSPOSE) {
91           P = Q_1d;
92           Q = Q_1d;
93         }
94         CeedInt           pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
95         const CeedScalar *interp_1d;
96 
97         CeedCallBackend(CeedBasisGetInterp1D(basis, &interp_1d));
98         if (impl->collo_grad_1d) {
99           CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)];
100           CeedScalar interp[num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)];
101 
102           // Interpolate to quadrature points (NoTranspose)
103           //  or Grad to quadrature points (Transpose)
104           for (CeedInt d = 0; d < dim; d++) {
105             CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? interp_1d : impl->collo_grad_1d), t_mode,
106                                                     (t_mode == CEED_TRANSPOSE) && (d > 0),
107                                                     (t_mode == CEED_NOTRANSPOSE ? (d == 0 ? u : tmp[d % 2]) : &u[d * num_qpts * num_comp * num_elem]),
108                                                     (t_mode == CEED_NOTRANSPOSE ? (d == dim - 1 ? interp : tmp[(d + 1) % 2]) : interp)));
109             pre /= P;
110             post *= Q;
111           }
112           // Grad to quadrature points (NoTranspose)
113           //  or Interpolate to nodes (Transpose)
114           P = Q_1d, Q = Q_1d;
115           if (t_mode == CEED_TRANSPOSE) {
116             P = Q_1d;
117             Q = P_1d;
118           }
119           pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
120           for (CeedInt d = 0; d < dim; d++) {
121             CeedCallBackend(CeedTensorContractApply(
122                 contract, pre, P, post, Q, (t_mode == CEED_NOTRANSPOSE ? impl->collo_grad_1d : interp_1d), t_mode,
123                 (t_mode == CEED_NOTRANSPOSE && apply_add) || (t_mode == CEED_TRANSPOSE && (d == dim - 1)),
124                 (t_mode == CEED_NOTRANSPOSE ? interp : (d == 0 ? interp : tmp[d % 2])),
125                 (t_mode == CEED_NOTRANSPOSE ? &v[d * num_qpts * num_comp * num_elem] : (d == dim - 1 ? v : tmp[(d + 1) % 2]))));
126             pre /= P;
127             post *= Q;
128           }
129         } else if (impl->has_collo_interp) {  // Qpts collocated with nodes
130           const CeedScalar *grad_1d;
131 
132           CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d));
133 
134           // Dim contractions, identity in other directions
135           CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
136 
137           for (CeedInt d = 0; d < dim; d++) {
138             CeedCallBackend(CeedTensorContractApply(contract, pre, P, post, Q, grad_1d, t_mode, add && (d > 0),
139                                                     t_mode == CEED_NOTRANSPOSE ? u : &u[d * num_comp * num_qpts * num_elem],
140                                                     t_mode == CEED_TRANSPOSE ? v : &v[d * num_comp * num_qpts * num_elem]));
141             pre /= P;
142             post *= Q;
143           }
144         } else {  // Underintegration, P > Q
145           const CeedScalar *grad_1d;
146 
147           CeedCallBackend(CeedBasisGetGrad1D(basis, &grad_1d));
148 
149           if (t_mode == CEED_TRANSPOSE) {
150             P = Q_1d;
151             Q = P_1d;
152           }
153           CeedScalar tmp[2][num_elem * num_comp * Q * CeedIntPow(P > Q ? P : Q, dim - 1)];
154 
155           // Dim**2 contractions, apply grad when pass == dim
156           for (CeedInt p = 0; p < dim; p++) {
157             CeedInt pre = num_comp * CeedIntPow(P, dim - 1), post = num_elem;
158 
159             for (CeedInt d = 0; d < dim; d++) {
160               CeedCallBackend(CeedTensorContractApply(
161                   contract, pre, P, post, Q, (p == d) ? grad_1d : interp_1d, t_mode, add && (d == dim - 1),
162                   (d == 0 ? (t_mode == CEED_NOTRANSPOSE ? u : &u[p * num_comp * num_qpts * num_elem]) : tmp[d % 2]),
163                   (d == dim - 1 ? (t_mode == CEED_TRANSPOSE ? v : &v[p * num_comp * num_qpts * num_elem]) : tmp[(d + 1) % 2])));
164               pre /= P;
165               post *= Q;
166             }
167           }
168         }
169       } break;
170       // Retrieve interpolation weights
171       case CEED_EVAL_WEIGHT: {
172         CeedInt           Q = Q_1d;
173         const CeedScalar *q_weight_1d;
174 
175         CeedCheck(t_mode == CEED_NOTRANSPOSE, ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
176         CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight_1d));
177         for (CeedInt d = 0; d < dim; d++) {
178           CeedInt pre = CeedIntPow(Q, dim - d - 1), post = CeedIntPow(Q, d);
179 
180           for (CeedInt i = 0; i < pre; i++) {
181             for (CeedInt j = 0; j < Q; j++) {
182               for (CeedInt k = 0; k < post; k++) {
183                 const CeedScalar w = q_weight_1d[j] * (d == 0 ? 1 : v[((i * Q + j) * post + k) * num_elem]);
184 
185                 for (CeedInt e = 0; e < num_elem; e++) v[((i * Q + j) * post + k) * num_elem + e] = w;
186               }
187             }
188           }
189         }
190       } break;
191       // LCOV_EXCL_START
192       case CEED_EVAL_DIV:
193       case CEED_EVAL_CURL:
194         return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
195       case CEED_EVAL_NONE:
196         return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
197         // LCOV_EXCL_STOP
198     }
199   } else {
200     // Non-tensor basis
201     CeedInt P = num_nodes, Q = num_qpts;
202 
203     switch (eval_mode) {
204       // Interpolate to/from quadrature points
205       case CEED_EVAL_INTERP: {
206         const CeedScalar *interp;
207 
208         CeedCallBackend(CeedBasisGetInterp(basis, &interp));
209         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, interp, t_mode, add, u, v));
210       } break;
211       // Evaluate the gradient to/from quadrature points
212       case CEED_EVAL_GRAD: {
213         const CeedScalar *grad;
214 
215         CeedCallBackend(CeedBasisGetGrad(basis, &grad));
216         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, grad, t_mode, add, u, v));
217       } break;
218       // Evaluate the divergence to/from the quadrature points
219       case CEED_EVAL_DIV: {
220         const CeedScalar *div;
221 
222         CeedCallBackend(CeedBasisGetDiv(basis, &div));
223         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, div, t_mode, add, u, v));
224       } break;
225       // Evaluate the curl to/from the quadrature points
226       case CEED_EVAL_CURL: {
227         const CeedScalar *curl;
228 
229         CeedCallBackend(CeedBasisGetCurl(basis, &curl));
230         CeedCallBackend(CeedTensorContractStridedApply(contract, num_comp, P, num_elem, q_comp, Q, curl, t_mode, add, u, v));
231       } break;
232       // Retrieve interpolation weights
233       case CEED_EVAL_WEIGHT: {
234         const CeedScalar *q_weight;
235 
236         CeedCheck(t_mode == CEED_NOTRANSPOSE, ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE");
237         CeedCallBackend(CeedBasisGetQWeights(basis, &q_weight));
238         for (CeedInt i = 0; i < num_qpts; i++) {
239           for (CeedInt e = 0; e < num_elem; e++) v[i * num_elem + e] = q_weight[i];
240         }
241       } break;
242       // LCOV_EXCL_START
243       case CEED_EVAL_NONE:
244         return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_NONE does not make sense in this context");
245         // LCOV_EXCL_STOP
246     }
247   }
248   if (U != CEED_VECTOR_NONE) {
249     CeedCallBackend(CeedVectorRestoreArrayRead(U, &u));
250   }
251   CeedCallBackend(CeedVectorRestoreArray(V, &v));
252   return CEED_ERROR_SUCCESS;
253 }
254 
255 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
256   CeedCallBackend(CeedBasisApplyCore_Ref(basis, false, num_elem, t_mode, eval_mode, U, V));
257   return CEED_ERROR_SUCCESS;
258 }
259 
260 static int CeedBasisApplyAdd_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
261   CeedCallBackend(CeedBasisApplyCore_Ref(basis, true, num_elem, t_mode, eval_mode, U, V));
262   return CEED_ERROR_SUCCESS;
263 }
264 
265 //------------------------------------------------------------------------------
266 // Basis Destroy Tensor
267 //------------------------------------------------------------------------------
268 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) {
269   CeedBasis_Ref *impl;
270 
271   CeedCallBackend(CeedBasisGetData(basis, &impl));
272   CeedCallBackend(CeedFree(&impl->collo_grad_1d));
273   CeedCallBackend(CeedFree(&impl));
274   return CEED_ERROR_SUCCESS;
275 }
276 
277 //------------------------------------------------------------------------------
278 // Basis Create Tensor
279 //------------------------------------------------------------------------------
280 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
281                                 const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis) {
282   Ceed               ceed, ceed_parent;
283   CeedBasis_Ref     *impl;
284   CeedTensorContract contract;
285 
286   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
287   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
288 
289   CeedCallBackend(CeedCalloc(1, &impl));
290   // Check for collocated interp
291   if (Q_1d == P_1d) {
292     bool has_collocated = true;
293 
294     for (CeedInt i = 0; i < P_1d; i++) {
295       has_collocated = has_collocated && (fabs(interp_1d[i + P_1d * i] - 1.0) < 1e-14);
296       for (CeedInt j = 0; j < P_1d; j++) {
297         if (j != i) has_collocated = has_collocated && (fabs(interp_1d[j + P_1d * i]) < 1e-14);
298       }
299     }
300     impl->has_collo_interp = has_collocated;
301   }
302   // Calculate collocated grad
303   if (Q_1d >= P_1d && !impl->has_collo_interp) {
304     CeedCallBackend(CeedMalloc(Q_1d * Q_1d, &impl->collo_grad_1d));
305     CeedCallBackend(CeedBasisGetCollocatedGrad(basis, impl->collo_grad_1d));
306   }
307   CeedCallBackend(CeedBasisSetData(basis, impl));
308 
309   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
310   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
311 
312   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
313   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
314   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroyTensor_Ref));
315   return CEED_ERROR_SUCCESS;
316 }
317 
318 //------------------------------------------------------------------------------
319 // Basis Create Non-Tensor H^1
320 //------------------------------------------------------------------------------
321 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *grad,
322                           const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
323   Ceed               ceed, ceed_parent;
324   CeedTensorContract contract;
325 
326   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
327   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
328 
329   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
330   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
331 
332   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
333   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
334   return CEED_ERROR_SUCCESS;
335 }
336 
337 //------------------------------------------------------------------------------
338 // Basis Create Non-Tensor H(div)
339 //------------------------------------------------------------------------------
340 int CeedBasisCreateHdiv_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div,
341                             const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
342   Ceed               ceed, ceed_parent;
343   CeedTensorContract contract;
344 
345   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
346   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
347 
348   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
349   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
350 
351   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
352   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
353   return CEED_ERROR_SUCCESS;
354 }
355 
356 //------------------------------------------------------------------------------
357 // Basis Create Non-Tensor H(curl)
358 //------------------------------------------------------------------------------
359 int CeedBasisCreateHcurl_Ref(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
360                              const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
361   Ceed               ceed, ceed_parent;
362   CeedTensorContract contract;
363 
364   CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
365   CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
366 
367   CeedCallBackend(CeedTensorContractCreate(ceed_parent, &contract));
368   CeedCallBackend(CeedBasisSetTensorContract(basis, contract));
369 
370   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApply_Ref));
371   CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAdd_Ref));
372   return CEED_ERROR_SUCCESS;
373 }
374 
375 //------------------------------------------------------------------------------
376