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