xref: /libCEED/backends/ref/ceed-ref-operator.c (revision 57c64913cb6f06c2f6a40825dd942ceba0d2fe12) 
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, &evecs[i+starti]);
61       CeedChk(ierr);
62     }
63 
64     switch(emode) {
65     case CEED_EVAL_NONE:
66       break; // No action
67     case CEED_EVAL_INTERP:
68       ncomp = qfields[i].ncomp;
69       ierr = CeedMalloc(Q*ncomp, &qdata_alloc[iq]); CeedChk(ierr);
70       qdata[i + starti] = qdata_alloc[iq];
71       iq++;
72       break;
73     case CEED_EVAL_GRAD:
74       ncomp = qfields[i].ncomp;
75       dim = ofields[i].basis->dim;
76       ierr = CeedMalloc(Q*ncomp*dim, &qdata_alloc[iq]); CeedChk(ierr);
77       qdata[i + starti] = qdata_alloc[iq];
78       iq++;
79       break;
80     case CEED_EVAL_WEIGHT: // Only on input fields
81       ierr = CeedMalloc(Q, &qdata_alloc[iq]); CeedChk(ierr);
82       ierr = CeedBasisApply(ofields[iq].basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT,
83                             NULL, qdata_alloc[iq]); CeedChk(ierr);
84       qdata[i] = qdata_alloc[iq];
85       indata[i] = qdata[i];
86       iq++;
87       break;
88     case CEED_EVAL_DIV:
89       break; // Not implimented
90     case CEED_EVAL_CURL:
91       break; // Not implimented
92     }
93   }
94   return 0;
95 }
96 
97 /*
98   CeedOperator needs to connect all the named fields (be they active or passive)
99   to the named inputs and outputs of its CeedQFunction.
100  */
101 static int CeedOperatorSetup_Ref(CeedOperator op) {
102   if (op->setupdone) return 0;
103   CeedOperator_Ref *impl = op->data;
104   CeedQFunction qf = op->qf;
105   CeedInt Q = op->numqpoints;
106   int ierr;
107 
108   // Count infield and outfield array sizes and evectors
109   impl->numein = qf->numinputfields;
110   for (CeedInt i=0; i<qf->numinputfields; i++) {
111     CeedEvalMode emode = qf->inputfields[i].emode;
112     impl->numqin += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD) + !!
113                      (emode & CEED_EVAL_WEIGHT);
114   }
115   impl->numeout = qf->numoutputfields;
116   for (CeedInt i=0; i<qf->numoutputfields; i++) {
117     CeedEvalMode emode = qf->outputfields[i].emode;
118     impl->numqout += !!(emode & CEED_EVAL_INTERP) + !!(emode & CEED_EVAL_GRAD);
119   }
120 
121   // Allocate
122   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->evecs); CeedChk(ierr);
123   ierr = CeedCalloc(impl->numein + impl->numeout, &impl->edata);
124   CeedChk(ierr);
125 
126   ierr = CeedCalloc(impl->numqin + impl->numqout, &impl->qdata_alloc);
127   CeedChk(ierr);
128   ierr = CeedCalloc(qf->numinputfields + qf->numoutputfields, &impl->qdata);
129   CeedChk(ierr);
130 
131   ierr = CeedCalloc(16, &impl->indata); CeedChk(ierr);
132   ierr = CeedCalloc(16, &impl->outdata); CeedChk(ierr);
133 
134   // Set up infield and outfield pointer arrays
135   // Infields
136   ierr = CeedOperatorSetupFields_Ref(qf->inputfields, op->inputfields,
137                                      impl->evecs, impl->qdata, impl->qdata_alloc,
138                                      impl->indata, 0, 0,
139                                      qf->numinputfields, Q); CeedChk(ierr);
140 
141   // Outfields
142   ierr = CeedOperatorSetupFields_Ref(qf->outputfields, op->outputfields,
143                                      impl->evecs, impl->qdata, impl->qdata_alloc,
144                                      impl->indata, qf->numinputfields,
145                                      impl->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       impl->outdata[i] =  impl->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 *impl = op->data;
163   CeedInt Q = op->numqpoints, elemsize;
164   int ierr;
165   CeedQFunction qf = op->qf;
166   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
167 
168   // Setup
169   ierr = CeedOperatorSetup_Ref(op); CeedChk(ierr);
170 
171   // Input Evecs and Restriction
172   for (CeedInt i=0; i<qf->numinputfields; i++) {
173     CeedEvalMode emode = qf->inputfields[i].emode;
174     if (emode == CEED_EVAL_WEIGHT) { // Skip
175     } else {
176       // Zero evec
177       ierr = CeedVectorSetValue(impl->evecs[i], 0.0); CeedChk(ierr);
178       // Active
179       if (op->inputfields[i].vec == CEED_VECTOR_ACTIVE) {
180         // Restrict
181         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
182                                         lmode, invec, impl->evecs[i],
183                                         request); CeedChk(ierr);
184         // Get evec
185         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
186                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
187       } else {
188         // Passive
189         // Restrict
190         ierr = CeedElemRestrictionApply(op->inputfields[i].Erestrict, CEED_NOTRANSPOSE,
191                                         lmode, op->inputfields[i].vec, impl->evecs[i],
192                                         request); CeedChk(ierr);
193         // Get evec
194         ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_HOST,
195                                       (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
196       }
197     }
198   }
199 
200   // Output Evecs
201   for (CeedInt i=0; i<qf->numoutputfields; i++) {
202     ierr = CeedVectorGetArray(impl->evecs[i+impl->numein], CEED_MEM_HOST,
203                               &impl->edata[i + qf->numinputfields]); CeedChk(ierr);
204   }
205 
206   // Loop through elements
207   for (CeedInt e=0; e<op->numelements; e++) {
208     // Input basis apply if needed
209     for (CeedInt i=0; i<qf->numinputfields; i++) {
210       // Get elemsize, emode, ncomp
211       elemsize = op->inputfields[i].Erestrict->elemsize;
212       CeedEvalMode emode = qf->inputfields[i].emode;
213       CeedInt ncomp = qf->inputfields[i].ncomp;
214       // Basis action
215       switch(emode) {
216       case CEED_EVAL_NONE:
217         impl->indata[i] = &impl->edata[i][e*Q*ncomp];
218         break;
219       case CEED_EVAL_INTERP:
220         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
221                               CEED_EVAL_INTERP, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
222         CeedChk(ierr);
223         impl->indata[i] = impl->qdata[i];
224         break;
225       case CEED_EVAL_GRAD:
226         ierr = CeedBasisApply(op->inputfields[i].basis, 1, CEED_NOTRANSPOSE,
227                               CEED_EVAL_GRAD, &impl->edata[i][e*elemsize*ncomp], impl->qdata[i]);
228         CeedChk(ierr);
229         impl->indata[i] = impl->qdata[i];
230         break;
231       case CEED_EVAL_WEIGHT:
232         break;  // No action
233       case CEED_EVAL_DIV:
234         break; // Not implimented
235       case CEED_EVAL_CURL:
236         break; // Not implimented
237       }
238     }
239     // Output pointers
240     for (CeedInt i=0; i<qf->numoutputfields; i++) {
241       CeedEvalMode emode = qf->outputfields[i].emode;
242       if (emode == CEED_EVAL_NONE) {
243         CeedInt ncomp = qf->outputfields[i].ncomp;
244         impl->outdata[i] = &impl->edata[i + qf->numinputfields][e*Q*ncomp];
245       }
246     }
247     // Q function
248     ierr = CeedQFunctionApply(op->qf, Q, (const CeedScalar * const*) impl->indata,
249                               impl->outdata); CeedChk(ierr);
250 
251     // Output basis apply if needed
252     for (CeedInt i=0; i<qf->numoutputfields; i++) {
253       // Get elemsize, emode, ncomp
254       elemsize = op->outputfields[i].Erestrict->elemsize;
255       CeedInt ncomp = qf->outputfields[i].ncomp;
256       CeedEvalMode emode = qf->outputfields[i].emode;
257       // Basis action
258       switch(emode) {
259       case CEED_EVAL_NONE:
260         break; // No action
261       case CEED_EVAL_INTERP:
262         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
263                               CEED_EVAL_INTERP, impl->outdata[i],
264                               &impl->edata[i + qf->numinputfields][e*elemsize*ncomp]); CeedChk(ierr);
265         break;
266       case CEED_EVAL_GRAD:
267         ierr = CeedBasisApply(op->outputfields[i].basis, 1, CEED_TRANSPOSE,
268                               CEED_EVAL_GRAD,
269                               impl->outdata[i], &impl->edata[i + qf->numinputfields][e*elemsize*ncomp]);
270         CeedChk(ierr);
271         break;
272       case CEED_EVAL_WEIGHT:
273         break; // Should not occur
274       case CEED_EVAL_DIV:
275         break; // Not implimented
276       case CEED_EVAL_CURL:
277         break; // Not implimented
278       }
279     }
280   }
281 
282   // Output restriction
283   for (CeedInt i=0; i<qf->numoutputfields; i++) {
284     // Restore evec
285     ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein],
286                                     &impl->edata[i + qf->numinputfields]); CeedChk(ierr);
287     // Active
288     if (op->outputfields[i].vec == CEED_VECTOR_ACTIVE) {
289       // Zero lvec
290       ierr = CeedVectorSetValue(outvec, 0.0); CeedChk(ierr);
291       // Restrict
292       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
293                                       lmode, impl->evecs[i+impl->numein], outvec, request); CeedChk(ierr);
294     } else {
295       // Passive
296       // Zero lvec
297       ierr = CeedVectorSetValue(op->outputfields[i].vec, 0.0); CeedChk(ierr);
298       // Restrict
299       ierr = CeedElemRestrictionApply(op->outputfields[i].Erestrict, CEED_TRANSPOSE,
300                                       lmode, impl->evecs[i+impl->numein], op->outputfields[i].vec,
301                                       request); CeedChk(ierr);
302     }
303   }
304 
305   // Restore input arrays
306   for (CeedInt i=0; i<qf->numinputfields; i++) {
307     CeedEvalMode emode = qf->inputfields[i].emode;
308     if (emode == CEED_EVAL_WEIGHT) { // Skip
309     } else {
310       ierr = CeedVectorRestoreArrayRead(impl->evecs[i],
311                                         (const CeedScalar **) &impl->edata[i]); CeedChk(ierr);
312     }
313   }
314 
315   return 0;
316 }
317 
318 int CeedOperatorCreate_Ref(CeedOperator op) {
319   CeedOperator_Ref *impl;
320   int ierr;
321 
322   ierr = CeedCalloc(1, &impl); CeedChk(ierr);
323   op->data = impl;
324   op->Destroy = CeedOperatorDestroy_Ref;
325   op->Apply = CeedOperatorApply_Ref;
326   return 0;
327 }
328