xref: /honee/src/qdata.c (revision 49c2c2da4bdb2e3c2e8a681550774f7954eb4237) 
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   *q_data_size = 0;
107   switch (dim) {
108     case 2:
109       switch (num_comp_x) {
110         case 2:
111           *q_data_size = 5;
112           break;
113         case 3:
114           *q_data_size = 7;
115           break;
116       }
117       break;
118     case 3:
119       *q_data_size = 10;
120       break;
121   }
122   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
123              "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
124   PetscFunctionReturn(PETSC_SUCCESS);
125 }
126 
127 /**
128  * @brief Create quadrature data for domain
129  *
130  * @param[in]  ceed          Ceed object quadrature data will be used with
131  * @param[in]  dm            DM where quadrature data would be used
132  * @param[in]  domain_label  DMLabel that quadrature data would be used one
133  * @param[in]  label_value   Value of label
134  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
135  * @param[in]  basis_x       CeedBasis of the coordinates
136  * @param[in]  x_coord       CeedVector of the coordinates
137  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
138  * @param[out] q_data        CeedVector of the quadrature data
139  * @param[out] q_data_size   Number of components of quadrature data
140  */
141 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
142                         CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
143   CeedQFunction       qf_setup = NULL;
144   CeedOperator        op_setup;
145   CeedVector          q_data_created;
146   CeedElemRestriction elem_restr_qd_created;
147   CeedInt             num_comp_x;
148   PetscInt            dim, height = 0;
149 
150   PetscFunctionBeginUser;
151   PetscCall(QDataGetNumComponents(dm, q_data_size));
152   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
153   PetscCall(DMGetDimension(dm, &dim));
154   switch (dim) {
155     case 2:
156       switch (num_comp_x) {
157         case 2:
158           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2d, Setup2d_loc, &qf_setup));
159           break;
160         case 3:
161           PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2D_3Dcoords, Setup2D_3Dcoords_loc, &qf_setup));
162           break;
163       }
164       break;
165     case 3:
166       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup, Setup_loc, &qf_setup));
167       break;
168   }
169   PetscCheck(qf_setup, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
170 
171   // -- Create QFunction for quadrature data
172   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup, 0));
173   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
174   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT));
175   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "surface qdata", *q_data_size, CEED_EVAL_NONE));
176 
177   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created));
178   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
179 
180   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
181   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
182   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
183   PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
184 
185   PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
186 
187   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
188 
189   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
190   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
191   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
192   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
193   PetscFunctionReturn(PETSC_SUCCESS);
194 }
195 
196 /**
197  * @brief Get number of components of quadrature data for boundary of domain
198  *
199  * @param[in]  dm          DM where quadrature data would be used
200  * @param[out] q_data_size Number of components of quadrature data
201  */
202 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) {
203   PetscInt dim;
204 
205   PetscFunctionBeginUser;
206   PetscCall(DMGetDimension(dm, &dim));
207   *q_data_size = 0;
208   switch (dim) {
209     case 2:
210       *q_data_size = 3;
211       break;
212     case 3:
213       *q_data_size = 10;
214       break;
215   }
216   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
217   PetscFunctionReturn(PETSC_SUCCESS);
218 }
219 
220 /**
221  * @brief Create quadrature data for boundary of domain
222  *
223  * @param[in]  ceed          Ceed object quadrature data will be used with
224  * @param[in]  dm            DM where quadrature data would be used
225  * @param[in]  domain_label  DMLabel that quadrature data would be used one
226  * @param[in]  label_value   Value of label
227  * @param[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
228  * @param[in]  basis_x       CeedBasis of the coordinates
229  * @param[in]  x_coord       CeedVector of the coordinates
230  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
231  * @param[out] q_data        CeedVector of the quadrature data
232  * @param[out] q_data_size   Number of components of quadrature data
233  */
234 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
235                                 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
236   CeedQFunction       qf_setup_sur = NULL;
237   CeedOperator        op_setup_sur;
238   CeedVector          q_data_created;
239   CeedElemRestriction elem_restr_qd_created;
240   CeedInt             num_comp_x;
241   PetscInt            dim, height = 1;
242 
243   PetscFunctionBeginUser;
244   PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
245   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
246   PetscCall(DMGetDimension(dm, &dim));
247   switch (dim) {
248     case 2:
249       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur));
250       break;
251     case 3:
252       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur));
253       break;
254   }
255   PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
256 
257   // -- Create QFunction for quadrature data
258   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0));
259   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
260   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
261   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE));
262 
263   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created));
264   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
265 
266   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur));
267   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
268   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
269   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
270 
271   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
272 
273   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
274 
275   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
276   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
277   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
278   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
279   PetscFunctionReturn(PETSC_SUCCESS);
280 }
281 
282 /**
283  * @brief Get number of components of quadrature data for boundary of domain
284  *
285  * @param[in]  dm          DM where quadrature data would be used
286  * @param[out] q_data_size Number of components of quadrature data
287  */
288 PetscErrorCode QDataBoundaryGradientGetNumComponents(DM dm, CeedInt *q_data_size) {
289   PetscInt dim;
290 
291   PetscFunctionBeginUser;
292   PetscCall(DMGetDimension(dm, &dim));
293   *q_data_size = 0;
294   switch (dim) {
295     case 2:
296       *q_data_size = 7;
297       break;
298     case 3:
299       *q_data_size = 13;
300       break;
301   }
302   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
303   PetscFunctionReturn(PETSC_SUCCESS);
304 }
305 
306 /**
307   @brief Compute `CeedOperator` surface gradient QData
308 
309   Collective across MPI processes.
310 
311   @param[in]  ceed          `Ceed` object
312   @param[in]  dm            `DMPlex` grid
313   @param[in]  domain_label  `DMLabel` for surface
314   @param[in]  label_value   `DMPlex` label value for surface
315   @param[in]  x_coord       `CeedVector` for coordinates
316   @param[out] elem_restr_qd `CeedElemRestriction` for QData
317   @param[out] q_data        `CeedVector` holding QData
318   @param[out] q_data_size   Number of QData components per quadrature point
319 **/
320 PetscErrorCode QDataBoundaryGradientGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedVector x_coord,
321                                         CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
322   PetscInt            dim;
323   const PetscInt      height_cell = 0, height_face = 1;
324   CeedInt             num_comp_x;
325   CeedElemRestriction elem_restr_x_cell, elem_restr_x_face, elem_restr_qd_created;
326   CeedVector          q_data_created;
327   CeedBasis           basis_x_cell_to_face, basis_x_face;
328   CeedQFunction       qf_setup_sur = NULL;
329   CeedOperator        op_setup_sur;
330 
331   PetscFunctionBeginUser;
332   PetscCall(CreateCoordinateBasisFromPlex(ceed, dm, domain_label, label_value, height_face, &basis_x_face));
333   PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, domain_label, label_value, label_value, &basis_x_cell_to_face));
334   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_cell, &elem_restr_x_cell));
335   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &elem_restr_x_face));
336 
337   PetscCall(QDataBoundaryGradientGetNumComponents(dm, q_data_size));
338   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height_face, *q_data_size, &elem_restr_qd_created));
339   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
340   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x_face, &num_comp_x));
341 
342   PetscCall(DMGetDimension(dm, &dim));
343   switch (dim) {
344     case 2:
345       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2DBoundaryGradient, Setup2DBoundaryGradient_loc, &qf_setup_sur));
346       break;
347     case 3:
348       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundaryGradient, SetupBoundaryGradient_loc, &qf_setup_sur));
349       break;
350   }
351   PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
352 
353   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX cell", num_comp_x * (dim - height_cell), CEED_EVAL_GRAD));
354   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX face", num_comp_x * (dim - height_face), CEED_EVAL_GRAD));
355   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
356   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "q data", *q_data_size, CEED_EVAL_NONE));
357 
358   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_sur));
359   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX cell", elem_restr_x_cell, basis_x_cell_to_face, CEED_VECTOR_ACTIVE));
360   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX face", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE));
361   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE));
362   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "q data", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
363 
364   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
365 
366   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
367 
368   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
369   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
370   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
371   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
372   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
373   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
374   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
375   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
376   PetscFunctionReturn(PETSC_SUCCESS);
377 }
378