xref: /honee/src/qdata.c (revision e77831d27f21b8ea2cfbf833b1f1fb3c0911dd7c) 
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 
4 #include <navierstokes.h>
5 
6 #include <petscsection.h>
7 #include "../qfunctions/setupgeo.h"
8 #include "../qfunctions/setupgeo2d.h"
9 
10 static CeedVector          *q_data_vecs;
11 static CeedElemRestriction *q_data_restrictions;
12 static PetscInt             num_q_data_stored;
13 
14 /**
15   @brief Get stored QData objects that match created objects, if present
16 
17   If created objects are not present, they are added to the storage and returned in the output
18 
19   Not Collective
20 
21   @param[in]  q_data_created        Vector with created QData
22   @param[in]  elem_restr_qd_created Restriction for created QData
23   @param[out] q_data_stored         Vector from storage matching QData
24   @param[out] elem_restr_qd_stored  Restriction from storage matching QData
25 **/
26 static PetscErrorCode QDataGetStored(CeedVector q_data_created, CeedElemRestriction elem_restr_qd_created, CeedVector *q_data_stored,
27                                      CeedElemRestriction *elem_restr_qd_stored) {
28   Ceed     ceed = CeedVectorReturnCeed(q_data_created);
29   CeedSize created_length, stored_length;
30   PetscInt q_data_stored_index = -1;
31 
32   PetscFunctionBeginUser;
33   PetscCallCeed(ceed, CeedVectorGetLength(q_data_created, &created_length));
34   for (PetscInt i = 0; i < num_q_data_stored; i++) {
35     CeedVector difference_cvec;
36     CeedScalar max_difference;
37 
38     PetscCallCeed(ceed, CeedVectorGetLength(q_data_vecs[0], &stored_length));
39     if (created_length != stored_length) continue;
40     PetscCallCeed(ceed, CeedVectorCreate(ceed, stored_length, &difference_cvec));
41     PetscCallCeed(ceed, CeedVectorCopy(q_data_vecs[i], difference_cvec));
42     PetscCallCeed(ceed, CeedVectorAXPY(difference_cvec, -1, q_data_created));
43     PetscCallCeed(ceed, CeedVectorNorm(difference_cvec, CEED_NORM_MAX, &max_difference));
44     PetscCallCeed(ceed, CeedVectorDestroy(&difference_cvec));
45     if (max_difference < 100 * CEED_EPSILON) {
46       q_data_stored_index = i;
47       break;
48     }
49   }
50 
51   if (q_data_stored_index == -1) {
52     q_data_stored_index = num_q_data_stored++;
53 
54     PetscCall(PetscRealloc(num_q_data_stored * sizeof(CeedVector), &q_data_vecs));
55     PetscCall(PetscRealloc(num_q_data_stored * sizeof(CeedElemRestriction), &q_data_restrictions));
56     q_data_vecs[q_data_stored_index]         = NULL;  // Must set to NULL for ReferenceCopy
57     q_data_restrictions[q_data_stored_index] = NULL;  // Must set to NULL for ReferenceCopy
58     PetscCallCeed(ceed, CeedVectorReferenceCopy(q_data_created, &q_data_vecs[q_data_stored_index]));
59     PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(elem_restr_qd_created, &q_data_restrictions[q_data_stored_index]));
60   }
61   *q_data_stored        = NULL;  // Must set to NULL for ReferenceCopy
62   *elem_restr_qd_stored = NULL;  // Must set to NULL for ReferenceCopy
63   PetscCallCeed(ceed, CeedVectorReferenceCopy(q_data_vecs[q_data_stored_index], q_data_stored));
64   PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(q_data_restrictions[q_data_stored_index], elem_restr_qd_stored));
65   PetscFunctionReturn(PETSC_SUCCESS);
66 }
67 
68 /**
69   @brief Clear stored QData objects
70 **/
71 PetscErrorCode QDataClearStoredData() {
72   PetscFunctionBeginUser;
73   for (PetscInt i = 0; i < num_q_data_stored; i++) {
74     Ceed ceed = CeedVectorReturnCeed(q_data_vecs[i]);
75 
76     PetscCallCeed(ceed, CeedVectorDestroy(&q_data_vecs[i]));
77     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&q_data_restrictions[i]));
78   }
79   PetscCall(PetscFree(q_data_vecs));
80   PetscCall(PetscFree(q_data_restrictions));
81   q_data_vecs         = NULL;
82   q_data_restrictions = NULL;
83   PetscFunctionReturn(PETSC_SUCCESS);
84 }
85 
86 /**
87  * @brief Get number of components of quadrature data for domain
88  *
89  * @param[in]  dm          DM where quadrature data would be used
90  * @param[out] q_data_size Number of components of quadrature data
91  */
92 PetscErrorCode QDataGetNumComponents(DM dm, CeedInt *q_data_size) {
93   PetscInt num_comp_x, dim;
94 
95   PetscFunctionBeginUser;
96   PetscCall(DMGetDimension(dm, &dim));
97   {  // Get number of coordinate components
98     DM           dm_coord;
99     PetscSection section_coord;
100     PetscInt     field = 0;  // Default field has the coordinates
101 
102     PetscCall(DMGetCoordinateDM(dm, &dm_coord));
103     PetscCall(DMGetLocalSection(dm_coord, &section_coord));
104     PetscCall(PetscSectionGetFieldComponents(section_coord, field, &num_comp_x));
105   }
106   switch (dim) {
107     case 2:
108       switch (num_comp_x) {
109         case 2:
110           *q_data_size = 5;
111           break;
112         case 3:
113           *q_data_size = 7;
114           break;
115         default:
116           SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
117                   "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
118           break;
119       }
120       break;
121     case 3:
122       *q_data_size = 10;
123       break;
124     default:
125       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
126               "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
127       break;
128   }
129   PetscFunctionReturn(PETSC_SUCCESS);
130 }
131 
132 /**
133  * @brief Create quadrature data for domain
134  *
135  * @param[in]  ceed          Ceed object quadrature data will be used with
136  * @param[in]  dm            DM where quadrature data would be used
137  * @param[in]  domain_label  DMLabel that quadrature data would be used one
138  * @param[in]  label_value   Value of label
139  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
140  * @param[in]  basis_x       CeedBasis of the coordinates
141  * @param[in]  x_coord       CeedVector of the coordinates
142  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
143  * @param[out] q_data        CeedVector of the quadrature data
144  * @param[out] q_data_size   Number of components of quadrature data
145  */
146 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
147                         CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
148   CeedQFunction       qf_setup = NULL;
149   CeedOperator        op_setup;
150   CeedVector          q_data_created;
151   CeedElemRestriction elem_restr_qd_created;
152   CeedInt             num_comp_x;
153   PetscInt            dim, height = 0;
154 
155   PetscFunctionBeginUser;
156   PetscCall(QDataGetNumComponents(dm, q_data_size));
157   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
158   PetscCall(DMGetDimension(dm, &dim));
159   switch (dim) {
160     case 2:
161       switch (num_comp_x) {
162         case 2:
163           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2d, Setup2d_loc, &qf_setup));
164           break;
165         case 3:
166           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2D_3Dcoords, Setup2D_3Dcoords_loc, &qf_setup));
167           break;
168       }
169       break;
170     case 3:
171       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup, Setup_loc, &qf_setup));
172       break;
173   }
174 
175   // -- Create QFunction for quadrature data
176   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup, 0));
177   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
178   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT));
179   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "surface qdata", *q_data_size, CEED_EVAL_NONE));
180 
181   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created));
182   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
183 
184   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
185   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
186   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
187   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
188 
189   PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
190 
191   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
192 
193   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
194   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
195   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
196   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
197   PetscFunctionReturn(PETSC_SUCCESS);
198 }
199 
200 /**
201  * @brief Get number of components of quadrature data for boundary of domain
202  *
203  * @param[in]  dm          DM where quadrature data would be used
204  * @param[out] q_data_size Number of components of quadrature data
205  */
206 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) {
207   PetscInt dim;
208 
209   PetscFunctionBeginUser;
210   PetscCall(DMGetDimension(dm, &dim));
211   switch (dim) {
212     case 2:
213       *q_data_size = 3;
214       break;
215     case 3:
216       *q_data_size = 10;
217       break;
218     default:
219       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "QDataBoundary not valid for DM of dimension %" PetscInt_FMT, dim);
220       break;
221   }
222   PetscFunctionReturn(PETSC_SUCCESS);
223 }
224 
225 /**
226  * @brief Create quadrature data for boundary of domain
227  *
228  * @param[in]  ceed          Ceed object quadrature data will be used with
229  * @param[in]  dm            DM where quadrature data would be used
230  * @param[in]  domain_label  DMLabel that quadrature data would be used one
231  * @param[in]  label_value   Value of label
232  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
233  * @param[in]  basis_x       CeedBasis of the coordinates
234  * @param[in]  x_coord       CeedVector of the coordinates
235  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
236  * @param[out] q_data        CeedVector of the quadrature data
237  * @param[out] q_data_size   Number of components of quadrature data
238  */
239 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
240                                 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
241   CeedQFunction       qf_setup_sur = NULL;
242   CeedOperator        op_setup_sur;
243   CeedVector          q_data_created;
244   CeedElemRestriction elem_restr_qd_created;
245   CeedInt             num_comp_x;
246   PetscInt            dim, height = 1;
247 
248   PetscFunctionBeginUser;
249   PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
250   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
251   PetscCall(DMGetDimension(dm, &dim));
252   switch (dim) {
253     case 2:
254       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur));
255       break;
256     case 3:
257       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur));
258       break;
259   }
260 
261   // -- Create QFunction for quadrature data
262   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0));
263   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
264   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
265   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE));
266 
267   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created));
268   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
269 
270   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur));
271   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
272   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
273   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
274 
275   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
276 
277   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
278 
279   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
280   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
281   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
282   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
283   PetscFunctionReturn(PETSC_SUCCESS);
284 }
285 
286 /**
287  * @brief Get number of components of quadrature data for boundary of domain
288  *
289  * @param[in]  dm          DM where quadrature data would be used
290  * @param[out] q_data_size Number of components of quadrature data
291  */
292 PetscErrorCode QDataBoundaryGradientGetNumComponents(DM dm, CeedInt *q_data_size) {
293   PetscInt dim;
294 
295   PetscFunctionBeginUser;
296   PetscCall(DMGetDimension(dm, &dim));
297   switch (dim) {
298     case 3:
299       *q_data_size = 13;
300       break;
301     default:
302       SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "QDataBoundary not valid for DM of dimension %" PetscInt_FMT, dim);
303       break;
304   }
305   PetscFunctionReturn(PETSC_SUCCESS);
306 }
307 
308 /**
309   @brief Compute `CeedOperator` surface gradient QData
310 
311   Collective across MPI processes.
312 
313   @param[in]  ceed          `Ceed` object
314   @param[in]  dm            `DMPlex` grid
315   @param[in]  domain_label  `DMLabel` for surface
316   @param[in]  label_value   `DMPlex` label value for surface
317   @param[in]  x_coord       `CeedVector` for coordinates
318   @param[out] elem_restr_qd `CeedElemRestriction` for QData
319   @param[out] q_data        `CeedVector` holding QData
320   @param[out] q_data_size   Number of QData components per quadrature point
321 **/
322 PetscErrorCode QDataBoundaryGradientGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedVector x_coord,
323                                         CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
324   PetscInt            dim;
325   const PetscInt      height_cell = 0, height_face = 1;
326   CeedInt             num_comp_x;
327   CeedElemRestriction elem_restr_x_cell, elem_restr_x_face, elem_restr_qd_created;
328   CeedVector          q_data_created;
329   CeedBasis           basis_x_cell_to_face, basis_x_face;
330   CeedQFunction       qf_setup_sur;
331   CeedOperator        op_setup_sur;
332 
333   PetscFunctionBeginUser;
334   PetscCall(CreateCoordinateBasisFromPlex(ceed, dm, domain_label, label_value, height_face, &basis_x_face));
335   PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, domain_label, label_value, label_value, &basis_x_cell_to_face));
336   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_cell, &elem_restr_x_cell));
337   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &elem_restr_x_face));
338 
339   PetscCall(QDataBoundaryGradientGetNumComponents(dm, q_data_size));
340   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height_face, *q_data_size, &elem_restr_qd_created));
341   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
342 
343   PetscCall(DMGetDimension(dm, &dim));
344   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x_face, &num_comp_x));
345 
346   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundaryGradient, SetupBoundaryGradient_loc, &qf_setup_sur));
347   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX cell", num_comp_x * (dim - height_cell), CEED_EVAL_GRAD));
348   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX face", num_comp_x * (dim - height_face), CEED_EVAL_GRAD));
349   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
350   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "q data", *q_data_size, CEED_EVAL_NONE));
351 
352   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_sur));
353   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX cell", elem_restr_x_cell, basis_x_cell_to_face, CEED_VECTOR_ACTIVE));
354   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX face", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE));
355   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE));
356   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "q data", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
357 
358   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
359 
360   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
361 
362   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
363   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
364   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
365   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
366   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
367   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
368   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
369   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
370   PetscFunctionReturn(PETSC_SUCCESS);
371 }
372