xref: /libCEED/backends/ref/ceed-ref-operator.c (revision a4b33194745f1db2c078785db78486bb2b14379e) 
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 starte,
52                                        CeedInt startq, CeedInt numfields, CeedInt Q) {
53   CeedInt dim, ierr, ie=starte, iq=startq, ncomp;
54 
55   // Loop over fields
56   for (CeedInt i=0; i<numfields; i++) {
57     if (ofields[i].Erestrict != CEED_RESTRICTION_IDENTITY) {
58       ierr = CeedElemRestrictionCreateVector(ofields[i].Erestrict, NULL, &evecs[ie]);
59       CeedChk(ierr);
60       ie++;
61     }
62     CeedEvalMode emode = qfields[i].emode;
63     switch(emode) {
64     case CEED_EVAL_NONE:
65       break; // No action
66     case CEED_EVAL_INTERP:
67       ncomp = qfields[i].ncomp;
68       ierr = CeedMalloc(Q*ncomp, &qdata_alloc[iq]); CeedChk(ierr);
69       qdata[i + starti] = qdata_alloc[iq];
70       iq++;
71       break;
72     case CEED_EVAL_GRAD:
73       ncomp = qfields[i].ncomp;
74       dim = ofields[i].basis->dim;
75       ierr = CeedMalloc(Q*ncomp*dim, &qdata_alloc[iq]); CeedChk(ierr);
76       qdata[i + starti] = qdata_alloc[iq];
77       iq++;
78       break;
79     case CEED_EVAL_WEIGHT: // Only on input fields
80       ierr = CeedMalloc(Q, &qdata_alloc[iq]); CeedChk(ierr);
81       ierr = CeedBasisApply(ofields[iq].basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
82                             NULL, qdata_alloc[iq]); CeedChk(ierr);
83       qdata[i] = qdata_alloc[iq];
84       indata[i] = qdata[i];
85       iq++;
86       break;
87     case CEED_EVAL_DIV:
88       break; // Not implimented
89     case CEED_EVAL_CURL:
90       break; // Not implimented
91     }
92   }
93   return 0;
94 }
95 
96 /*
97   CeedOperator needs to connect all the named fields (be they active or passive)
98   to the named inputs and outputs of its CeedQFunction.
99  */
100 static int CeedOperatorSetup_Ref(CeedOperator op) {
101   if (op->setupdone) return 0;
102   CeedOperator_Ref *opref = op->data;
103   CeedQFunction qf = op->qf;
104   CeedInt Q = op->numqpoints;
105   int ierr;
106 
107   // Count infield and outfield array sizes and evectors
108   for (CeedInt i=0; i<qf->numinputfields; i++) {
109     CeedEvalMode emode = qf->inputfields[i].emode;
110     opref->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) + !!
111                      (emode & CEED_EVAL_WEIGHT);
112     opref->numein +=
113       (op->inputfields[i].Erestrict !=
114        CEED_RESTRICTION_IDENTITY); // Need E-vector when non-identity restriction exists
115   }
116   for (CeedInt i=0; i<qf->numoutputfields; i++) {
117     CeedEvalMode emode = qf->outputfields[i].emode;
118     opref->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD);
119     opref->numeout += (op->outputfields[i].Erestrict != CEED_RESTRICTION_IDENTITY);
120   }
121 
122   // Allocate
123   ierr = CeedCalloc(opref->numein + opref->numeout, &opref->evecs); CeedChk(ierr);
124   ierr = CeedCalloc(qf->numinputfields + qf->numoutputfields, &opref->edata);
125   CeedChk(ierr);
126 
127   ierr = CeedCalloc(opref->numqin + opref->numqout, &opref->qdata_alloc);
128   CeedChk(ierr);
129   ierr = CeedCalloc(qf->numinputfields + qf->numoutputfields, &opref->qdata);
130   CeedChk(ierr);
131 
132   ierr = CeedCalloc(16, &opref->indata); CeedChk(ierr);
133   ierr = CeedCalloc(16, &opref->outdata); CeedChk(ierr);
134 
135   // Set up infield and outfield pointer arrays
136   // Infields
137   ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields,
138                                      opref->evecs, opref->qdata, opref->qdata_alloc,
139                                      opref->indata, 0, 0, 0,
140                                      qf->numinputfields, Q); CeedChk(ierr);
141 
142   // Outfields
143   ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields,
144                                      opref->evecs, opref->qdata, opref->qdata_alloc,
145                                      opref->indata, qf->numinputfields, opref->numein,
146                                      opref->numqin, qf->numoutputfields, Q); CeedChk(ierr);
147 
148   // Output Qvecs
149   for (CeedInt i=0; i<qf->numoutputfields; i++) {
150     CeedEvalMode emode = qf->outputfields[i].emode;
151     if (emode != CEED_EVAL_NONE) {
152       opref->outdata[i] =  opref->qdata[i + qf->numinputfields];
153     }
154   }
155 
156   op->setupdone = 1;
157 
158   return 0;
159 }
160 
161 static int CeedOperatorApply_Ref(CeedOperator op, CeedVector invec,
162                                  CeedVector outvec, CeedRequest *request) {
163   CeedOperator_Ref *opref = op->data;
164   CeedInt Q = op->numqpoints, elemsize;
165   int ierr;
166   CeedQFunction qf = op->qf;
167   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
168   CeedScalar *vec_temp;
169 
170   // Setup
171   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
172 
173   // Input Evecs and Restriction
174   for (CeedInt i=0,iein=0; i<qf->numinputfields; i++) {
175     // No Restriction
176     if (op->inputfields[i].Erestrict == CEED_RESTRICTION_IDENTITY) {
177       CeedEvalMode emode = qf->inputfields[i].emode;
178       if (emode & CEED_EVAL_WEIGHT) {
179       } else {
180         // Active
181         if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
182           ierr = CeedVectorGetArrayRead(invec, CEED_MEM_HOST,
183                                         (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
184           // Passive
185         } else {
186           ierr = CeedVectorGetArrayRead(op->inputfields[i].vec, CEED_MEM_HOST,
187                                         (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
188         }
189       }
190     } else {
191       // Restriction
192       // Zero evec
193       ierr = CeedVectorGetArray(opref->evecs[iein], CEED_MEM_HOST, &vec_temp);
194       CeedChk(ierr);
195       for (CeedInt j=0; j<opref->evecs[iein]->length; j++)
196         vec_temp[j] = 0.;
197       ierr = CeedVectorRestoreArray(opref->evecs[iein], &vec_temp); CeedChk(ierr);
198       // Active
199       if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
200         // Restrict
201         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
202                                         lmode, invec, opref->evecs[iein],
203                                         request); CeedChk(ierr);
204         // Get evec
205         ierr = CeedVectorGetArrayRead(opref->evecs[iein], CEED_MEM_HOST,
206                                       (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
207         iein++;
208       } else {
209         // Passive
210         // Restrict
211         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
212                                         lmode, op->inputfields[i].vec, opref->evecs[iein],
213                                         request); CeedChk(ierr);
214         // Get evec
215         ierr = CeedVectorGetArrayRead(opref->evecs[iein], CEED_MEM_HOST,
216                                       (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
217         iein++;
218       }
219     }
220   }
221 
222   // Output Evecs
223   for (CeedInt i=0,ieout=opref->numein; i<qf->numoutputfields; i++) {
224     // No Restriction
225     if (op->outputfields[i].Erestrict == CEED_RESTRICTION_IDENTITY) {
226       // Active
227       if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
228         ierr = CeedVectorGetArray(outvec, CEED_MEM_HOST,
229                                   &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
230       } else {
231         // Passive
232         ierr = CeedVectorGetArray(op->outputfields[i].vec, CEED_MEM_HOST,
233                                   &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
234       }
235     } else {
236       // Restriction
237       ierr = CeedVectorGetArray(opref->evecs[ieout], CEED_MEM_HOST,
238                                 &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
239       ieout++;
240     }
241   }
242 
243   // Loop through elements
244   for (CeedInt e=0; e<op->numelements; e++) {
245     // Input basis apply if needed
246     for (CeedInt i=0; i<qf->numinputfields; i++) {
247       // Get elemsize
248       if (op->inputfields[i].Erestrict != CEED_RESTRICTION_IDENTITY) {
249         elemsize = op->inputfields[i].Erestrict->elemsize;
250       } else {
251         elemsize = Q;
252       }
253       // Get emode, ncomp
254       CeedEvalMode emode = qf->inputfields[i].emode;
255       CeedInt ncomp = qf->inputfields[i].ncomp;
256       // Basis action
257       switch(emode) {
258       case CEED_EVAL_NONE:
259         opref->indata[i] = &opref->edata[i][e*Q*ncomp];
260         break;
261       case CEED_EVAL_INTERP:
262         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
263                               CEED_EVAL_INTERP, &opref->edata[i][e*elemsize*ncomp], opref->qdata[i]);
264         CeedChk(ierr);
265         opref->indata[i] = opref->qdata[i];
266         break;
267       case CEED_EVAL_GRAD:
268         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
269                               CEED_EVAL_GRAD, &opref->edata[i][e*elemsize*ncomp], opref->qdata[i]);
270         CeedChk(ierr);
271         opref->indata[i] = opref->qdata[i];
272         break;
273       case CEED_EVAL_WEIGHT:
274         break;  // No action
275       case CEED_EVAL_DIV:
276         break; // Not implimented
277       case CEED_EVAL_CURL:
278         break; // Not implimented
279       }
280     }
281     // Output pointers
282     for (CeedInt i=0; i<qf->numoutputfields; i++) {
283       CeedEvalMode emode = qf->outputfields[i].emode;
284       if (emode == CEED_EVAL_NONE) {
285         CeedInt ncomp = qf->outputfields[i].ncomp;
286         opref->outdata[i] = &opref->edata[i + qf->numinputfields][e*Q*ncomp];
287       }
288     }
289     // Q function
290     ierr = CeedQFunctionApply(op->qf, Q, (const CeedScalar * const*) opref->indata,
291                               opref->outdata); CeedChk(ierr);
292 
293     // Output basis apply if needed
294     for (CeedInt i=0; i<qf->numoutputfields; i++) {
295       // Get elemsize
296       if (op->outputfields[i].Erestrict != CEED_RESTRICTION_IDENTITY) {
297         elemsize = op->outputfields[i].Erestrict->elemsize;
298       } else {
299         elemsize = Q;
300       }
301       // Get emode, ncomp
302       CeedInt ncomp = qf->outputfields[i].ncomp;
303       CeedEvalMode emode = qf->outputfields[i].emode;
304       // Basis action
305       switch(emode) {
306       case CEED_EVAL_NONE:
307         break; // No action
308       case CEED_EVAL_INTERP:
309         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
310                               CEED_EVAL_INTERP, opref->outdata[i],
311                               &opref->edata[i + qf->numinputfields][e*elemsize*ncomp]); CeedChk(ierr);
312         break;
313       case CEED_EVAL_GRAD:
314         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
315                               CEED_EVAL_GRAD,
316                               opref->outdata[i], &opref->edata[i + qf->numinputfields][e*elemsize*ncomp]);
317         CeedChk(ierr);
318         break;
319       case CEED_EVAL_WEIGHT:
320         break; // Should not occur
321       case CEED_EVAL_DIV:
322         break; // Not implimented
323       case CEED_EVAL_CURL:
324         break; // Not implimented
325       }
326     }
327   }
328 
329   // Output restriction
330   for (CeedInt i=0,ieout=opref->numein; i<qf->numoutputfields; i++) {
331     // No Restriction
332     if (op->outputfields[i].Erestrict == CEED_RESTRICTION_IDENTITY) {
333       // Active
334       if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
335         ierr = CeedVectorRestoreArray(outvec, &opref->edata[i + qf->numinputfields]);
336         CeedChk(ierr);
337       } else {
338         // Passive
339         ierr = CeedVectorRestoreArray(op->outputfields[i].vec,
340                                       &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
341       }
342     } else {
343       // Restriction
344       // Active
345       if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
346         // Restore evec
347         ierr = CeedVectorRestoreArray(opref->evecs[ieout],
348                                       &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
349         // Zero lvec
350         ierr = CeedVectorGetArray(outvec, CEED_MEM_HOST, &vec_temp); CeedChk(ierr);
351         for (CeedInt j=0; j<outvec->length; j++)
352           vec_temp[j] = 0.;
353         ierr = CeedVectorRestoreArray(outvec, &vec_temp); CeedChk(ierr);
354         // Restrict
355         ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
356                                         lmode, opref->evecs[ieout], outvec, request); CeedChk(ierr);
357         ieout++;
358       } else {
359         // Passive
360         // Restore evec
361         ierr = CeedVectorRestoreArray(opref->evecs[ieout],
362                                       &opref->edata[i + qf->numinputfields]); CeedChk(ierr);
363         // Zero lvec
364         ierr = CeedVectorGetArray(op->outputfields[i].vec, CEED_MEM_HOST, &vec_temp);
365         CeedChk(ierr);
366         for (CeedInt j=0; j<op->outputfields[i].vec->length; j++)
367           vec_temp[j] = 0.;
368         ierr = CeedVectorRestoreArray(op->outputfields[i].vec, &vec_temp);
369         CeedChk(ierr);
370         // Restrict
371         ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
372                                         lmode, opref->evecs[ieout], op->outputfields[i].vec,
373                                         request); CeedChk(ierr);
374         ieout++;
375       }
376     }
377   }
378 
379   // Restore input arrays
380   for (CeedInt i=0,iein=0; i<qf->numinputfields; i++) {
381     // No Restriction
382     if (op->inputfields[i].Erestrict == CEED_RESTRICTION_IDENTITY) {
383       CeedEvalMode emode = qf->inputfields[i].emode;
384       if (emode & CEED_EVAL_WEIGHT) {
385       } else {
386         // Active
387         if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
388           ierr = CeedVectorRestoreArrayRead(invec,
389                                             (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
390           // Passive
391         } else {
392           ierr = CeedVectorRestoreArrayRead(op->inputfields[i].vec,
393                                             (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
394         }
395       }
396     } else {
397       // Restriction
398       ierr = CeedVectorRestoreArrayRead(opref->evecs[iein],
399                                         (const CeedScalar **) &opref->edata[i]); CeedChk(ierr);
400       iein++;
401     }
402   }
403 
404   return 0;
405 }
406 
407 int CeedOperatorCreate_Ref(CeedOperator op) {
408   CeedOperator_Ref *impl;
409   int ierr;
410 
411   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
412   op->data = impl;
413   op->Destroy = CeedOperatorDestroy_Ref;
414   op->Apply = CeedOperatorApply_Ref;
415   return 0;
416 }
417