xref: /libCEED/backends/ref/ceed-ref-operator.c (revision d863ab9ba6a0d47c58445a35d35b36efba07fc93) 
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->numqin+impl->numqout; i++) {
32     ierr = CeedFree(&impl->qdata_alloc[i]); CeedChk(ierr);
33   }
34   ierr = CeedFree(&impl->qdata_alloc); CeedChk(ierr);
35   ierr = CeedFree(&impl->qdata); CeedChk(ierr);
36 
37   ierr = CeedFree(&impl->indata); CeedChk(ierr);
38   ierr = CeedFree(&impl->outdata); CeedChk(ierr);
39 
40   ierr = CeedFree(&op->data); CeedChk(ierr);
41   return 0;
42 }
43 
44 /*
45   Setup infields or outfields
46  */
47 static int CeedOperatorSetupFields_Ref(struct CeedQFunctionField qfields[16],
48                                        struct CeedOperatorField ofields[16],
49                                        CeedVector *evecs, CeedScalar **qdata,
50                                        CeedScalar **qdata_alloc, CeedScalar **indata,
51                                        CeedInt starti, CeedInt startq,
52                                        CeedInt numfields, CeedInt Q) {
53   CeedInt dim, ierr, iq=startq, ncomp;
54 
55   // Loop over fields
56   for (CeedInt i=0; i<numfields; i++) {
57     CeedEvalMode emode = qfields[i].emode;
58 
59     if (emode != CEED_EVAL_WEIGHT) {
60       ierr = CeedElemRestrictionCreateVector(ofields[i].Erestrict, NULL,
61                                              &evecs[i+starti]);
62       CeedChk(ierr);
63     }
64 
65     switch(emode) {
66     case CEED_EVAL_NONE:
67       break; // No action
68     case CEED_EVAL_INTERP:
69       ncomp = qfields[i].ncomp;
70       ierr = CeedMalloc(Q*ncomp, &qdata_alloc[iq]); CeedChk(ierr);
71       qdata[i + starti] = qdata_alloc[iq];
72       iq++;
73       break;
74     case CEED_EVAL_GRAD:
75       ncomp = qfields[i].ncomp;
76       ierr = CeedBasisGetDimension(ofields[i].basis, &dim); CeedChk(ierr);
77       ierr = CeedMalloc(Q*ncomp*dim, &qdata_alloc[iq]); CeedChk(ierr);
78       qdata[i + starti] = qdata_alloc[iq];
79       iq++;
80       break;
81     case CEED_EVAL_WEIGHT: // Only on input fields
82       ierr = CeedMalloc(Q, &qdata_alloc[iq]); CeedChk(ierr);
83       ierr = CeedBasisApply(ofields[iq].basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
84                             NULL, qdata_alloc[iq]); CeedChk(ierr);
85       qdata[i] = qdata_alloc[iq];
86       indata[i] = qdata[i];
87       iq++;
88       break;
89     case CEED_EVAL_DIV:
90       break; // Not implimented
91     case CEED_EVAL_CURL:
92       break; // Not implimented
93     }
94   }
95   return 0;
96 }
97 
98 /*
99   CeedOperator needs to connect all the named fields (be they active or passive)
100   to the named inputs and outputs of its CeedQFunction.
101  */
102 static int CeedOperatorSetup_Ref(CeedOperator op) {
103   int ierr;
104   bool setupdone;
105   ierr = CeedOperatorGetSetupStatus(op, &setupdone); CeedChk(ierr);
106   if (setupdone) return 0;
107   CeedOperator_Ref *impl;
108   ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr);
109   CeedQFunction qf;
110   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
111   CeedInt Q, numinputfields, numoutputfields;
112   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
113   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
114   CeedChk(ierr);
115 
116   // Count infield and outfield array sizes and evectors
117   impl->numein = numinputfields;
118   for (CeedInt i=0; i<numinputfields; i++) {
119     CeedEvalMode emode = qf->inputfields[i].emode;
120     impl->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) +
121                     !!(emode & CEED_EVAL_WEIGHT);
122   }
123   impl->numeout = numoutputfields;
124   for (CeedInt i=0; i<numoutputfields; i++) {
125     CeedEvalMode emode = qf->outputfields[i].emode;
126     impl->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD);
127   }
128 
129   // Allocate
130   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->evecs); CeedChk(ierr);
131   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->edata);
132   CeedChk(ierr);
133 
134   ierr = CeedCalloc(impl->numqin + impl->numqout, &impl->qdata_alloc);
135   CeedChk(ierr);
136   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->qdata);
137   CeedChk(ierr);
138 
139   ierr = CeedCalloc(16, &impl->indata); CeedChk(ierr);
140   ierr = CeedCalloc(16, &impl->outdata); CeedChk(ierr);
141 
142   // Set up infield and outfield pointer arrays
143   // Infields
144   ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields,
145                                      impl->evecs, impl->qdata, impl->qdata_alloc,
146                                      impl->indata, 0, 0,
147                                      numinputfields, Q); CeedChk(ierr);
148 
149   // Outfields
150   ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields,
151                                      impl->evecs, impl->qdata, impl->qdata_alloc,
152                                      impl->indata, numinputfields,
153                                      impl->numqin, numoutputfields, Q); CeedChk(ierr);
154 
155   // Input Qvecs
156   for (CeedInt i=0; i<numinputfields; i++) {
157     CeedEvalMode emode = qf->inputfields[i].emode;
158     if ((emode != CEED_EVAL_NONE) && (emode != CEED_EVAL_WEIGHT))
159       impl->indata[i] =  impl->qdata[i];
160   }
161   // Output Qvecs
162   for (CeedInt i=0; i<numoutputfields; i++) {
163     CeedEvalMode emode = qf->outputfields[i].emode;
164     if (emode != CEED_EVAL_NONE)
165       impl->outdata[i] =  impl->qdata[i + numinputfields];
166   }
167 
168   ierr = CeedOperatorSetSetupDone(op); CeedChk(ierr);
169 
170   return 0;
171 }
172 
173 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
174                                  CeedVector outvec, CeedRequest *request) {
175   int ierr;
176   CeedOperator_Ref *impl;
177   ierr = CeedOperatorGetData(op, (void*)&impl); CeedChk(ierr);
178   CeedQFunction qf;
179   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
180   CeedInt Q, numelements, elemsize, numinputfields, numoutputfields;
181   ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChk(ierr);
182   ierr = CeedOperatorGetNumElements(op, &numelements); CeedChk(ierr);
183   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
184   CeedChk(ierr);
185   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
186 
187   // Setup
188   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
189 
190   // Input Evecs and Restriction
191   for (CeedInt i=0; i<numinputfields; i++) {
192     CeedEvalMode emode = qf->inputfields[i].emode;
193     if (emode == CEED_EVAL_WEIGHT) { // Skip
194     } else {
195       // Active
196       if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
197         // Restrict
198         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
199                                         lmode, invec, impl->evecs[i],
200                                         request); CeedChk(ierr);
201         // Get evec
202         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
203                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
204       } else {
205         // Passive
206         // Restrict
207         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
208                                         lmode, op->inputfields[i].vec, impl->evecs[i],
209                                         request); CeedChk(ierr);
210         // Get evec
211         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
212                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
213       }
214     }
215   }
216 
217   // Output Evecs
218   for (CeedInt i=0; i<numoutputfields; i++) {
219     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
220                               &impl->edata[i + numinputfields]); CeedChk(ierr);
221   }
222 
223   // Loop through elements
224   for (CeedInt e=0; e<numelements; e++) {
225     // Input basis apply if needed
226     for (CeedInt i=0; i<numinputfields; i++) {
227       // Get elemsize, emode, ncomp
228       ierr = CeedElemRestrictionGetElementSize(
229                op->inputfields[i].Erestrict, &elemsize); CeedChk(ierr);
230       CeedEvalMode emode = qf->inputfields[i].emode;
231       CeedInt ncomp = qf->inputfields[i].ncomp;
232       // Basis action
233       switch(emode) {
234       case CEED_EVAL_NONE:
235         impl->indata[i] = &impl->edata[i][e*Q*ncomp];
236         break;
237       case CEED_EVAL_INTERP:
238         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
239                               CEED_EVAL_INTERP, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
240         CeedChk(ierr);
241         break;
242       case CEED_EVAL_GRAD:
243         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
244                               CEED_EVAL_GRAD, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
245         CeedChk(ierr);
246         break;
247       case CEED_EVAL_WEIGHT:
248         break;  // No action
249       case CEED_EVAL_DIV:
250         break; // Not implimented
251       case CEED_EVAL_CURL:
252         break; // Not implimented
253       }
254     }
255     // Output pointers
256     for (CeedInt i=0; i<numoutputfields; i++) {
257       CeedEvalMode emode = qf->outputfields[i].emode;
258       if (emode == CEED_EVAL_NONE) {
259         CeedInt ncomp = qf->outputfields[i].ncomp;
260         impl->outdata[i] = &impl->edata[i + numinputfields][e*Q*ncomp];
261       }
262     }
263     // Q function
264     ierr = CeedQFunctionApply(qf, Q, (const CeedScalar * const*) impl->indata,
265                               impl->outdata); CeedChk(ierr);
266 
267     // Output basis apply if needed
268     for (CeedInt i=0; i<numoutputfields; i++) {
269       // Get elemsize, emode, ncomp
270       ierr = CeedElemRestrictionGetElementSize(
271                op->outputfields[i].Erestrict, &elemsize); CeedChk(ierr);
272       CeedInt ncomp = qf->outputfields[i].ncomp;
273       CeedEvalMode emode = qf->outputfields[i].emode;
274       // Basis action
275       switch(emode) {
276       case CEED_EVAL_NONE:
277         break; // No action
278       case CEED_EVAL_INTERP:
279         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
280                               CEED_EVAL_INTERP, impl->outdata[i],
281                               &impl->edata[i + numinputfields][e*elemsize*ncomp]);
282         CeedChk(ierr);
283         break;
284       case CEED_EVAL_GRAD:
285         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
286                               CEED_EVAL_GRAD,
287                               impl->outdata[i], &impl->edata[i + numinputfields][e*elemsize*ncomp]);
288         CeedChk(ierr);
289         break;
290       case CEED_EVAL_WEIGHT: {
291         Ceed ceed;
292         ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
293         return CeedError(ceed, 1,
294                          "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
295         break; // Should not occur
296       }
297       case CEED_EVAL_DIV:
298         break; // Not implimented
299       case CEED_EVAL_CURL:
300         break; // Not implimented
301       }
302     }
303   }
304 
305   // Zero lvecs
306   ierr = CeedVectorSetValue(outvec, 0.0); CeedChk(ierr);
307   for (CeedInt i=0; i<numoutputfields; i++)
308     if (op->outputfields[i].vec != CEED_VECTOR_ACTIVE) {
309       ierr = CeedVectorSetValue(op->outputfields[i].vec, 0.0); CeedChk(ierr);
310     }
311 
312   // Output restriction
313   for (CeedInt i=0; i<numoutputfields; i++) {
314     // Restore evec
315     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
316                                   &impl->edata[i + numinputfields]); CeedChk(ierr);
317     // Active
318     if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
319       // Restrict
320       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
321                                       lmode, impl->evecs[i+impl->numein], outvec, request); CeedChk(ierr);
322     } else {
323       // Passive
324       // Restrict
325       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
326                                       lmode, impl->evecs[i+impl->numein], op->outputfields[i].vec,
327                                       request); CeedChk(ierr);
328     }
329   }
330 
331   // Restore input arrays
332   for (CeedInt i=0; i<numinputfields; i++) {
333     CeedEvalMode emode = qf->inputfields[i].emode;
334     if (emode == CEED_EVAL_WEIGHT) { // Skip
335     } else {
336       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
337                                         (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
338     }
339   }
340 
341   return 0;
342 }
343 
344 int CeedOperatorCreate_Ref(CeedOperator op) {
345   int ierr;
346   CeedOperator_Ref *impl;
347 
348   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
349   op->data = impl;
350   op->Destroy = CeedOperatorDestroy_Ref;
351   op->Apply = CeedOperatorApply_Ref;
352   return 0;
353 }
354