xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 783c99b3f03c04aed3140873345e25fc3bab2ee4) 
1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
3 // All Rights reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 #include <ceed-impl.h>
18 #include <string.h>
19 #include "ceed-ref.h"
20 
21 static int CeedOperatorDestroy_Ref(CeedOperator op) {
22   CeedOperator_Ref *impl = op->data;
23   int 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       dim = ofields[i].basis->dim;
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   if (op->setupdone) return 0;
104   CeedOperator_Ref *impl = op->data;
105   CeedQFunction qf = op->qf;
106   CeedInt Q = op->numqpoints, numinputfields, numoutputfields;
107   int ierr;
108 
109   // Count infield and outfield array sizes and evectors
110   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); CeedChk(ierr);
111   impl->numein = numinputfields;
112   for (CeedInt i=0; i<numinputfields; i++) {
113     CeedEvalMode emode = qf->inputfields[i].emode;
114     impl->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) +
115                     !!(emode & CEED_EVAL_WEIGHT);
116   }
117   impl->numeout = numoutputfields;
118   for (CeedInt i=0; i<numoutputfields; i++) {
119     CeedEvalMode emode = qf->outputfields[i].emode;
120     impl->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD);
121   }
122 
123   // Allocate
124   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->evecs); CeedChk(ierr);
125   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->edata);
126   CeedChk(ierr);
127 
128   ierr = CeedCalloc(impl->numqin + impl->numqout, &impl->qdata_alloc);
129   CeedChk(ierr);
130   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->qdata);
131   CeedChk(ierr);
132 
133   ierr = CeedCalloc(16, &impl->indata); CeedChk(ierr);
134   ierr = CeedCalloc(16, &impl->outdata); CeedChk(ierr);
135 
136   // Set up infield and outfield pointer arrays
137   // Infields
138   ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields,
139                                      impl->evecs, impl->qdata, impl->qdata_alloc,
140                                      impl->indata, 0, 0,
141                                      numinputfields, Q); CeedChk(ierr);
142 
143   // Outfields
144   ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields,
145                                      impl->evecs, impl->qdata, impl->qdata_alloc,
146                                      impl->indata, numinputfields,
147                                      impl->numqin, numoutputfields, Q); CeedChk(ierr);
148 
149   // Input Qvecs
150   for (CeedInt i=0; i<numinputfields; i++) {
151     CeedEvalMode emode = qf->inputfields[i].emode;
152     if ((emode != CEED_EVAL_NONE) && (emode != CEED_EVAL_WEIGHT))
153       impl->indata[i] =  impl->qdata[i];
154   }
155   // Output Qvecs
156   for (CeedInt i=0; i<numoutputfields; i++) {
157     CeedEvalMode emode = qf->outputfields[i].emode;
158     if (emode != CEED_EVAL_NONE)
159       impl->outdata[i] =  impl->qdata[i + numinputfields];
160   }
161 
162   op->setupdone = 1;
163 
164   return 0;
165 }
166 
167 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
168                                  CeedVector outvec, CeedRequest *request) {
169   CeedOperator_Ref *impl = op->data;
170   CeedInt Q = op->numqpoints, elemsize, numinputfields, numoutputfields;
171   int ierr;
172   CeedQFunction qf = op->qf;
173   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
174 
175   // Setup
176   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
177   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); CeedChk(ierr);
178 
179   // Input Evecs and Restriction
180   for (CeedInt i=0; i<numinputfields; i++) {
181     CeedEvalMode emode = qf->inputfields[i].emode;
182     if (emode == CEED_EVAL_WEIGHT) { // Skip
183     } else {
184       // Active
185       if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
186         // Restrict
187         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
188                                         lmode, invec, impl->evecs[i],
189                                         request); CeedChk(ierr);
190         // Get evec
191         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
192                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
193       } else {
194         // Passive
195         // Restrict
196         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
197                                         lmode, op->inputfields[i].vec, impl->evecs[i],
198                                         request); CeedChk(ierr);
199         // Get evec
200         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
201                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
202       }
203     }
204   }
205 
206   // Output Evecs
207   for (CeedInt i=0; i<numoutputfields; i++) {
208     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
209                               &impl->edata[i + numinputfields]); CeedChk(ierr);
210   }
211 
212   // Loop through elements
213   for (CeedInt e=0; e<op->numelements; e++) {
214     // Input basis apply if needed
215     for (CeedInt i=0; i<numinputfields; i++) {
216       // Get elemsize, emode, ncomp
217       elemsize = op->inputfields[i].Erestrict->elemsize;
218       CeedEvalMode emode = qf->inputfields[i].emode;
219       CeedInt ncomp = qf->inputfields[i].ncomp;
220       // Basis action
221       switch(emode) {
222       case CEED_EVAL_NONE:
223         impl->indata[i] = &impl->edata[i][e*Q*ncomp];
224         break;
225       case CEED_EVAL_INTERP:
226         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
227                               CEED_EVAL_INTERP, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
228         CeedChk(ierr);
229         break;
230       case CEED_EVAL_GRAD:
231         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
232                               CEED_EVAL_GRAD, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
233         CeedChk(ierr);
234         break;
235       case CEED_EVAL_WEIGHT:
236         break;  // No action
237       case CEED_EVAL_DIV:
238         break; // Not implimented
239       case CEED_EVAL_CURL:
240         break; // Not implimented
241       }
242     }
243     // Output pointers
244     for (CeedInt i=0; i<numoutputfields; i++) {
245       CeedEvalMode emode = qf->outputfields[i].emode;
246       if (emode == CEED_EVAL_NONE) {
247         CeedInt ncomp = qf->outputfields[i].ncomp;
248         impl->outdata[i] = &impl->edata[i + numinputfields][e*Q*ncomp];
249       }
250     }
251     // Q function
252     ierr = CeedQFunctionApply(op->qf, Q, (const CeedScalar * const*) impl->indata,
253                               impl->outdata); CeedChk(ierr);
254 
255     // Output basis apply if needed
256     for (CeedInt i=0; i<numoutputfields; i++) {
257       // Get elemsize, emode, ncomp
258       elemsize = op->outputfields[i].Erestrict->elemsize;
259       CeedInt ncomp = qf->outputfields[i].ncomp;
260       CeedEvalMode emode = qf->outputfields[i].emode;
261       // Basis action
262       switch(emode) {
263       case CEED_EVAL_NONE:
264         break; // No action
265       case CEED_EVAL_INTERP:
266         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
267                               CEED_EVAL_INTERP, impl->outdata[i],
268                               &impl->edata[i + numinputfields][e*elemsize*ncomp]);
269         CeedChk(ierr);
270         break;
271       case CEED_EVAL_GRAD:
272         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
273                               CEED_EVAL_GRAD,
274                               impl->outdata[i], &impl->edata[i + numinputfields][e*elemsize*ncomp]);
275         CeedChk(ierr);
276         break;
277       case CEED_EVAL_WEIGHT:
278         return CeedError(op->ceed, 1,
279                          "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
280         break; // Should not occur
281       case CEED_EVAL_DIV:
282         break; // Not implimented
283       case CEED_EVAL_CURL:
284         break; // Not implimented
285       }
286     }
287   }
288 
289   // Zero lvecs
290   ierr = CeedVectorSetValue(outvec, 0.0); CeedChk(ierr);
291   for (CeedInt i=0; i<qf->numoutputfields; i++)
292     if (op->outputfields[i].vec != CEED_VECTOR_ACTIVE) {
293       ierr = CeedVectorSetValue(op->outputfields[i].vec, 0.0); CeedChk(ierr);
294     }
295 
296   // Output restriction
297   for (CeedInt i=0; i<numoutputfields; i++) {
298     // Restore evec
299     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
300                                   &impl->edata[i + numinputfields]); CeedChk(ierr);
301     // Active
302     if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
303       // Restrict
304       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
305                                       lmode, impl->evecs[i+impl->numein], outvec, request); CeedChk(ierr);
306     } else {
307       // Passive
308       // Restrict
309       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
310                                       lmode, impl->evecs[i+impl->numein], op->outputfields[i].vec,
311                                       request); CeedChk(ierr);
312     }
313   }
314 
315   // Restore input arrays
316   for (CeedInt i=0; i<numinputfields; i++) {
317     CeedEvalMode emode = qf->inputfields[i].emode;
318     if (emode == CEED_EVAL_WEIGHT) { // Skip
319     } else {
320       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
321                                         (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
322     }
323   }
324 
325   return 0;
326 }
327 
328 int CeedOperatorCreate_Ref(CeedOperator op) {
329   CeedOperator_Ref *impl;
330   int ierr;
331 
332   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
333   op->data = impl;
334   op->Destroy = CeedOperatorDestroy_Ref;
335   op->Apply = CeedOperatorApply_Ref;
336   return 0;
337 }
338