xref: /libCEED/backends/ref/ceed-ref-operator.c (revision aedaa0e5fd3a3e03ad33ad8a6308ac527f4f900e) 
1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3 // All Rights reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 #include <string.h>
18 #include "ceed-ref.h"
19 
20 static int CeedOperatorDestroy_Ref(CeedOperator op) {
21   int ierr;
22   CeedOperator_Ref *impl;
23   ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr);
24 
25   for (CeedInt i=0; i<impl->numein+impl->numeout; i++) {
26     ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChk(ierr);
27   }
28   ierr = CeedFree(&impl->evecs); CeedChk(ierr);
29   ierr = CeedFree(&impl->edata); CeedChk(ierr);
30 
31   for (CeedInt i=0; i<impl->numein; i++) {
32     ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChk(ierr);
33   }
34   ierr = CeedFree(&impl->qvecsin); CeedChk(ierr);
35 
36   for (CeedInt i=0; i<impl->numeout; i++) {
37     ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChk(ierr);
38   }
39   ierr = CeedFree(&impl->qvecsout); CeedChk(ierr);
40 
41   ierr = CeedFree(&impl); CeedChk(ierr);
42   return 0;
43 }
44 
45 /*
46   Setup infields or outfields
47  */
48 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op,
49                                        bool inOrOut, CeedVector *evecs,
50                                        CeedVector *qvecs, CeedInt starte,
51                                        CeedInt numfields, CeedInt Q) {
52   CeedInt dim, ierr, ncomp;
53   Ceed ceed;
54   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
55   CeedBasis basis;
56   CeedElemRestriction Erestrict;
57   CeedOperatorField *opfields;
58   CeedQFunctionField *qffields;
59   if (inOrOut) {
60     ierr = CeedOperatorGetFields(op, NULL, &opfields);
61     CeedChk(ierr);
62     ierr = CeedQFunctionGetFields(qf, NULL, &qffields);
63     CeedChk(ierr);
64   } else {
65     ierr = CeedOperatorGetFields(op, &opfields, NULL);
66     CeedChk(ierr);
67     ierr = CeedQFunctionGetFields(qf, &qffields, NULL);
68     CeedChk(ierr);
69   }
70 
71   // Loop over fields
72   for (CeedInt i=0; i<numfields; i++) {
73     CeedEvalMode emode;
74     ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChk(ierr);
75 
76     if (emode != CEED_EVAL_WEIGHT) {
77       ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &Erestrict);
78       CeedChk(ierr);
79       ierr = CeedElemRestrictionCreateVector(Erestrict, NULL,
80                                              &evecs[i+starte]);
81       CeedChk(ierr);
82     }
83 
84     switch(emode) {
85     case CEED_EVAL_NONE:
86       ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp);
87       CeedChk(ierr);
88       ierr = CeedVectorCreate(ceed, Q*ncomp, &qvecs[i]); CeedChk(ierr);
89       break;
90     case CEED_EVAL_INTERP:
91       ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp);
92       CeedChk(ierr);
93       ierr = CeedVectorCreate(ceed, Q*ncomp, &qvecs[i]); CeedChk(ierr);
94       break;
95     case CEED_EVAL_GRAD:
96       ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr);
97       ierr = CeedQFunctionFieldGetNumComponents(qffields[i], &ncomp);
98       ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr);
99       ierr = CeedVectorCreate(ceed, Q*ncomp*dim, &qvecs[i]); CeedChk(ierr);
100       break;
101     case CEED_EVAL_WEIGHT: // Only on input fields
102       ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChk(ierr);
103       ierr = CeedVectorCreate(ceed, Q, &qvecs[i]); CeedChk(ierr);
104       ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
105                             NULL, qvecs[i]); CeedChk(ierr);
106       break;
107     case CEED_EVAL_DIV:
108       break; // Not implimented
109     case CEED_EVAL_CURL:
110       break; // Not implimented
111     }
112   }
113   return 0;
114 }
115 
116 /*
117   CeedOperator needs to connect all the named fields (be they active or passive)
118   to the named inputs and outputs of its CeedQFunction.
119  */
120 static int CeedOperatorSetup_Ref(CeedOperator op) {
121   int ierr;
122   bool setupdone;
123   ierr = CeedOperatorGetSetupStatus(op, &setupdone); CeedChk(ierr);
124   if (setupdone) return 0;
125   Ceed ceed;
126   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
127   CeedOperator_Ref *impl;
128   ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr);
129   CeedQFunction qf;
130   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
131   CeedInt Q, numinputfields, numoutputfields;
132   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
133   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
134   CeedChk(ierr);
135   CeedOperatorField *opinputfields, *opoutputfields;
136   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
137   CeedChk(ierr);
138   CeedQFunctionField *qfinputfields, *qfoutputfields;
139   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
140   CeedChk(ierr);
141 
142   // Allocate
143   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs);
144   CeedChk(ierr);
145   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->edata);
146   CeedChk(ierr);
147 
148   ierr = CeedCalloc(16, &impl->qvecsin); CeedChk(ierr);
149   ierr = CeedCalloc(16, &impl->qvecsout); CeedChk(ierr);
150 
151   impl->numein = numinputfields; impl->numeout = numoutputfields;
152 
153   // Set up infield and outfield evecs and qvecs
154   // Infields
155   ierr = CeedOperatorSetupFields_Ref(qf, op, 0,
156                                      impl->evecs, impl->qvecsin, 0,
157                                      numinputfields, Q);
158   CeedChk(ierr);
159 
160   // Outfields
161   ierr = CeedOperatorSetupFields_Ref(qf, op, 1,
162                                      impl->evecs, impl->qvecsout,
163                                      numinputfields, numoutputfields, Q);
164   CeedChk(ierr);
165 
166   // Temporary Vector
167   ierr = CeedVectorCreate(ceed, 0, &impl->tempvec); CeedChk(ierr);
168 
169   ierr = CeedOperatorSetSetupDone(op); CeedChk(ierr);
170 
171   return 0;
172 }
173 
174 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
175                                  CeedVector outvec, CeedRequest *request) {
176   int ierr;
177   CeedOperator_Ref *impl;
178   ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr);
179   CeedQFunction qf;
180   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
181   CeedInt Q, numelements, elemsize, numinputfields, numoutputfields, ncomp;
182   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
183   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
184   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
185   CeedChk(ierr);
186   CeedTransposeMode lmode;
187   CeedOperatorField *opinputfields, *opoutputfields;
188   ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields);
189   CeedChk(ierr);
190   CeedQFunctionField *qfinputfields, *qfoutputfields;
191   ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields);
192   CeedChk(ierr);
193   CeedEvalMode emode;
194   CeedVector vec;
195   CeedBasis basis;
196   CeedElemRestriction Erestrict;
197 
198   // Setup
199   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
200 
201   // Input Evecs and Restriction
202   for (CeedInt i=0; i<numinputfields; i++) {
203     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
204     CeedChk(ierr);
205     if (emode == CEED_EVAL_WEIGHT) { // Skip
206     } else {
207       // Get input vector
208       ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChk(ierr);
209       if (vec == CEED_VECTOR_ACTIVE)
210         vec = invec;
211       // Restrict
212       ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
213       CeedChk(ierr);
214       ierr = CeedOperatorFieldGetLMode(opinputfields[i], &lmode); CeedChk(ierr);
215       ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE,
216                                       lmode, vec, impl->evecs[i],
217                                       request); CeedChk(ierr);
218       // Get evec
219       ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
220                                     (const CeedScalar **) &impl->edata[i]);
221       CeedChk(ierr);
222     }
223   }
224 
225   // Output Evecs
226   for (CeedInt i=0; i<numoutputfields; i++) {
227     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
228                               &impl->edata[i + numinputfields]); CeedChk(ierr);
229   }
230 
231   // Loop through elements
232   for (CeedInt e=0; e<numelements; e++) {
233     // Input basis apply if needed
234     for (CeedInt i=0; i<numinputfields; i++) {
235       // Get elemsize, emode, ncomp
236       ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict);
237       CeedChk(ierr);
238       ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
239       CeedChk(ierr);
240       ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
241       CeedChk(ierr);
242       ierr = CeedQFunctionFieldGetNumComponents(qfinputfields[i], &ncomp);
243       CeedChk(ierr);
244       // Basis action
245       switch(emode) {
246       case CEED_EVAL_NONE:
247         ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_HOST,
248                                   CEED_USE_POINTER,
249                                   &impl->edata[i][e*Q*ncomp]); CeedChk(ierr);
250         break;
251       case CEED_EVAL_INTERP:
252         ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
253         ierr = CeedVectorSetArray(impl->tempvec, CEED_MEM_HOST,
254                                   CEED_USE_POINTER,
255                                   &impl->edata[i][e*elemsize*ncomp]);
256         CeedChk(ierr);
257         ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
258                               CEED_EVAL_INTERP, impl->tempvec,
259                               impl->qvecsin[i]); CeedChk(ierr);
260         break;
261       case CEED_EVAL_GRAD:
262         ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChk(ierr);
263         ierr = CeedVectorSetArray(impl->tempvec, CEED_MEM_HOST,
264                                   CEED_USE_POINTER,
265                                   &impl->edata[i][e*elemsize*ncomp]);
266         CeedChk(ierr);
267         ierr = CeedBasisApply(basis, 1, CEED_NOTRANSPOSE,
268                               CEED_EVAL_GRAD, impl->tempvec,
269                               impl->qvecsin[i]); CeedChk(ierr);
270         break;
271       case CEED_EVAL_WEIGHT:
272         break;  // No action
273       case CEED_EVAL_DIV:
274         break; // Not implimented
275       case CEED_EVAL_CURL:
276         break; // Not implimented
277       }
278     }
279     // Output pointers
280     for (CeedInt i=0; i<numoutputfields; i++) {
281       ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
282       CeedChk(ierr);
283       if (emode == CEED_EVAL_NONE) {
284         ierr = CeedQFunctionFieldGetNumComponents(qfoutputfields[i], &ncomp);
285         CeedChk(ierr);
286         ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_HOST,
287                                   CEED_USE_POINTER,
288                                   &impl->edata[i + numinputfields][e*Q*ncomp]);
289         CeedChk(ierr);
290       }
291     }
292     // Q function
293     ierr = CeedQFunctionApply(qf, Q, impl->qvecsin, impl->qvecsout); CeedChk(ierr);
294 
295     // Output basis apply if needed
296     for (CeedInt i=0; i<numoutputfields; i++) {
297       // Get elemsize, emode, ncomp
298       ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
299       CeedChk(ierr);
300       ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize);
301       CeedChk(ierr);
302       ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode);
303       CeedChk(ierr);
304       ierr = CeedQFunctionFieldGetNumComponents(qfoutputfields[i], &ncomp);
305       CeedChk(ierr);
306       // Basis action
307       switch(emode) {
308       case CEED_EVAL_NONE:
309         break; // No action
310       case CEED_EVAL_INTERP:
311         ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
312         CeedChk(ierr);
313         ierr = CeedVectorSetArray(impl->tempvec, CEED_MEM_HOST,
314                                   CEED_USE_POINTER,
315                                   &impl->edata[i + numinputfields][e*elemsize*ncomp]);
316         ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
317                               CEED_EVAL_INTERP, impl->qvecsout[i],
318                               impl->tempvec); CeedChk(ierr);
319         break;
320       case CEED_EVAL_GRAD:
321         ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis);
322         CeedChk(ierr);
323         ierr = CeedVectorSetArray(impl->tempvec, CEED_MEM_HOST,
324                                   CEED_USE_POINTER,
325                                   &impl->edata[i + numinputfields][e*elemsize*ncomp]);
326         ierr = CeedBasisApply(basis, 1, CEED_TRANSPOSE,
327                               CEED_EVAL_GRAD, impl->qvecsout[i],
328                               impl->tempvec); CeedChk(ierr);
329         break;
330       case CEED_EVAL_WEIGHT: {
331         Ceed ceed;
332         ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
333         return CeedError(ceed, 1,
334                          "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
335         break; // Should not occur
336       }
337       case CEED_EVAL_DIV:
338         break; // Not implimented
339       case CEED_EVAL_CURL:
340         break; // Not implimented
341       }
342     }
343   }
344 
345   // Zero lvecs
346   for (CeedInt i=0; i<numoutputfields; i++) {
347     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
348     if (vec == CEED_VECTOR_ACTIVE)
349       vec = outvec;
350     ierr = CeedVectorSetValue(vec, 0.0); CeedChk(ierr);
351   }
352 
353   // Output restriction
354   for (CeedInt i=0; i<numoutputfields; i++) {
355     // Restore evec
356     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
357                                   &impl->edata[i + numinputfields]);
358     CeedChk(ierr);
359     // Get output vector
360     ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); CeedChk(ierr);
361     // Active
362     if (vec == CEED_VECTOR_ACTIVE)
363       vec = outvec;
364     // Restrict
365     ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict);
366     CeedChk(ierr);
367     ierr = CeedOperatorFieldGetLMode(opoutputfields[i], &lmode); CeedChk(ierr);
368     ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE,
369                                     lmode, impl->evecs[i+impl->numein], vec,
370                                     request); CeedChk(ierr);
371   }
372 
373   // Restore input arrays
374   for (CeedInt i=0; i<numinputfields; i++) {
375     ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode);
376     CeedChk(ierr);
377     if (emode == CEED_EVAL_WEIGHT) { // Skip
378     } else {
379       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
380                                         (const CeedScalar **) &impl->edata[i]);
381       CeedChk(ierr);
382     }
383   }
384 
385   return 0;
386 }
387 
388 int CeedOperatorCreate_Ref(CeedOperator op) {
389   int ierr;
390   Ceed ceed;
391   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
392   CeedOperator_Ref *impl;
393 
394   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
395   ierr = CeedOperatorSetData(op, (void*)&impl);
396 
397   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Apply",
398                                 CeedOperatorApply_Ref); CeedChk(ierr);
399   ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",
400                                 CeedOperatorDestroy_Ref); CeedChk(ierr);
401   return 0;
402 }
403