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