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