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