xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 8da8c1b2815f330b68fd86c6d12e19570da8b93e) 
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);
111   CeedChk(ierr);
112   impl->numein = numinputfields;
113   for (CeedInt i=0; i<numinputfields; i++) {
114     CeedEvalMode emode = qf->inputfields[i].emode;
115     impl->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) +
116                     !!(emode & CEED_EVAL_WEIGHT);
117   }
118   impl->numeout = numoutputfields;
119   for (CeedInt i=0; i<numoutputfields; i++) {
120     CeedEvalMode emode = qf->outputfields[i].emode;
121     impl->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD);
122   }
123 
124   // Allocate
125   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->evecs); CeedChk(ierr);
126   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->edata);
127   CeedChk(ierr);
128 
129   ierr = CeedCalloc(impl->numqin + impl->numqout, &impl->qdata_alloc);
130   CeedChk(ierr);
131   ierr = CeedCalloc(numinputfields + numoutputfields, &impl->qdata);
132   CeedChk(ierr);
133 
134   ierr = CeedCalloc(16, &impl->indata); CeedChk(ierr);
135   ierr = CeedCalloc(16, &impl->outdata); CeedChk(ierr);
136 
137   // Set up infield and outfield pointer arrays
138   // Infields
139   ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields,
140                                      impl->evecs, impl->qdata, impl->qdata_alloc,
141                                      impl->indata, 0, 0,
142                                      numinputfields, Q); CeedChk(ierr);
143 
144   // Outfields
145   ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields,
146                                      impl->evecs, impl->qdata, impl->qdata_alloc,
147                                      impl->indata, numinputfields,
148                                      impl->numqin, numoutputfields, Q); CeedChk(ierr);
149 
150   // Input Qvecs
151   for (CeedInt i=0; i<numinputfields; i++) {
152     CeedEvalMode emode = qf->inputfields[i].emode;
153     if ((emode != CEED_EVAL_NONE) && (emode != CEED_EVAL_WEIGHT))
154       impl->indata[i] =  impl->qdata[i];
155   }
156   // Output Qvecs
157   for (CeedInt i=0; i<numoutputfields; i++) {
158     CeedEvalMode emode = qf->outputfields[i].emode;
159     if (emode != CEED_EVAL_NONE)
160       impl->outdata[i] =  impl->qdata[i + numinputfields];
161   }
162 
163   op->setupdone = 1;
164 
165   return 0;
166 }
167 
168 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
169                                  CeedVector outvec, CeedRequest *request) {
170   CeedOperator_Ref *impl = op->data;
171   CeedInt Q = op->numqpoints, elemsize, numinputfields, numoutputfields;
172   int ierr;
173   CeedQFunction qf = op->qf;
174   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
175 
176   // Setup
177   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
178   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
179   CeedChk(ierr);
180 
181   // Input Evecs and Restriction
182   for (CeedInt i=0; i<numinputfields; i++) {
183     CeedEvalMode emode = qf->inputfields[i].emode;
184     if (emode == CEED_EVAL_WEIGHT) { // Skip
185     } else {
186       // Active
187       if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
188         // Restrict
189         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
190                                         lmode, invec, impl->evecs[i],
191                                         request); CeedChk(ierr);
192         // Get evec
193         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
194                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
195       } else {
196         // Passive
197         // Restrict
198         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
199                                         lmode, op->inputfields[i].vec, 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       }
205     }
206   }
207 
208   // Output Evecs
209   for (CeedInt i=0; i<numoutputfields; i++) {
210     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
211                               &impl->edata[i + numinputfields]); CeedChk(ierr);
212   }
213 
214   // Loop through elements
215   for (CeedInt e=0; e<op->numelements; e++) {
216     // Input basis apply if needed
217     for (CeedInt i=0; i<numinputfields; i++) {
218       // Get elemsize, emode, ncomp
219       elemsize = op->inputfields[i].Erestrict->elemsize;
220       CeedEvalMode emode = qf->inputfields[i].emode;
221       CeedInt ncomp = qf->inputfields[i].ncomp;
222       // Basis action
223       switch(emode) {
224       case CEED_EVAL_NONE:
225         impl->indata[i] = &impl->edata[i][e*Q*ncomp];
226         break;
227       case CEED_EVAL_INTERP:
228         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
229                               CEED_EVAL_INTERP, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
230         CeedChk(ierr);
231         break;
232       case CEED_EVAL_GRAD:
233         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
234                               CEED_EVAL_GRAD, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
235         CeedChk(ierr);
236         break;
237       case CEED_EVAL_WEIGHT:
238         break;  // No action
239       case CEED_EVAL_DIV:
240         break; // Not implimented
241       case CEED_EVAL_CURL:
242         break; // Not implimented
243       }
244     }
245     // Output pointers
246     for (CeedInt i=0; i<numoutputfields; i++) {
247       CeedEvalMode emode = qf->outputfields[i].emode;
248       if (emode == CEED_EVAL_NONE) {
249         CeedInt ncomp = qf->outputfields[i].ncomp;
250         impl->outdata[i] = &impl->edata[i + numinputfields][e*Q*ncomp];
251       }
252     }
253     // Q function
254     ierr = CeedQFunctionApply(op->qf, Q, (const CeedScalar * const*) impl->indata,
255                               impl->outdata); CeedChk(ierr);
256 
257     // Output basis apply if needed
258     for (CeedInt i=0; i<numoutputfields; i++) {
259       // Get elemsize, emode, ncomp
260       elemsize = op->outputfields[i].Erestrict->elemsize;
261       CeedInt ncomp = qf->outputfields[i].ncomp;
262       CeedEvalMode emode = qf->outputfields[i].emode;
263       // Basis action
264       switch(emode) {
265       case CEED_EVAL_NONE:
266         break; // No action
267       case CEED_EVAL_INTERP:
268         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
269                               CEED_EVAL_INTERP, impl->outdata[i],
270                               &impl->edata[i + numinputfields][e*elemsize*ncomp]);
271         CeedChk(ierr);
272         break;
273       case CEED_EVAL_GRAD:
274         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
275                               CEED_EVAL_GRAD,
276                               impl->outdata[i], &impl->edata[i + numinputfields][e*elemsize*ncomp]);
277         CeedChk(ierr);
278         break;
279       case CEED_EVAL_WEIGHT:
280         return CeedError(op->ceed, 1,
281                          "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
282         break; // Should not occur
283       case CEED_EVAL_DIV:
284         break; // Not implimented
285       case CEED_EVAL_CURL:
286         break; // Not implimented
287       }
288     }
289   }
290 
291   // Zero lvecs
292   ierr = CeedVectorSetValue(outvec, 0.0); CeedChk(ierr);
293   for (CeedInt i=0; i<qf->numoutputfields; i++)
294     if (op->outputfields[i].vec != CEED_VECTOR_ACTIVE) {
295       ierr = CeedVectorSetValue(op->outputfields[i].vec, 0.0); CeedChk(ierr);
296     }
297 
298   // Output restriction
299   for (CeedInt i=0; i<numoutputfields; i++) {
300     // Restore evec
301     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
302                                   &impl->edata[i + numinputfields]); CeedChk(ierr);
303     // Active
304     if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
305       // Restrict
306       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
307                                       lmode, impl->evecs[i+impl->numein], outvec, request); CeedChk(ierr);
308     } else {
309       // Passive
310       // Restrict
311       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
312                                       lmode, impl->evecs[i+impl->numein], op->outputfields[i].vec,
313                                       request); CeedChk(ierr);
314     }
315   }
316 
317   // Restore input arrays
318   for (CeedInt i=0; i<numinputfields; i++) {
319     CeedEvalMode emode = qf->inputfields[i].emode;
320     if (emode == CEED_EVAL_WEIGHT) { // Skip
321     } else {
322       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
323                                         (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
324     }
325   }
326 
327   return 0;
328 }
329 
330 int CeedOperatorCreate_Ref(CeedOperator op) {
331   CeedOperator_Ref *impl;
332   int ierr;
333 
334   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
335   op->data = impl;
336   op->Destroy = CeedOperatorDestroy_Ref;
337   op->Apply = CeedOperatorApply_Ref;
338   return 0;
339 }
340