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