xref: /honee/src/qdata.c (revision 07126c9a638dba5b0bb1e3c5c6d0702032cbfb5e) 
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 <dm-utils.h>
7 #include <petscsection.h>
8 #include "../qfunctions/setupgeo.h"
9 #include "../qfunctions/setupgeo2d.h"
10 
11 static CeedVector          *q_data_vecs;
12 static CeedElemRestriction *q_data_restrictions;
13 static PetscInt             num_q_data_stored;
14 
15 /**
16   @brief Get stored QData objects that match created objects, if present
17 
18   If created objects are not present, they are added to the storage and returned in the output
19 
20   Not Collective
21 
22   @param[in]  q_data_created        Vector with created QData
23   @param[in]  elem_restr_qd_created Restriction for created QData
24   @param[out] q_data_stored         Vector from storage matching QData
25   @param[out] elem_restr_qd_stored  Restriction from storage matching QData
26 **/
27 static PetscErrorCode QDataGetStored(CeedVector q_data_created, CeedElemRestriction elem_restr_qd_created, CeedVector *q_data_stored,
28                                      CeedElemRestriction *elem_restr_qd_stored) {
29   Ceed     ceed = CeedVectorReturnCeed(q_data_created);
30   CeedSize created_length, stored_length;
31   PetscInt q_data_stored_index = -1;
32 
33   PetscFunctionBeginUser;
34   PetscCallCeed(ceed, CeedVectorGetLength(q_data_created, &created_length));
35   for (PetscInt i = 0; i < num_q_data_stored; i++) {
36     CeedVector difference_cvec;
37     CeedScalar max_difference;
38 
39     PetscCallCeed(ceed, CeedVectorGetLength(q_data_vecs[i], &stored_length));
40     if (created_length != stored_length) continue;
41     PetscCallCeed(ceed, CeedVectorCreate(ceed, stored_length, &difference_cvec));
42     PetscCallCeed(ceed, CeedVectorCopy(q_data_vecs[i], difference_cvec));
43     PetscCallCeed(ceed, CeedVectorAXPY(difference_cvec, -1, q_data_created));
44     PetscCallCeed(ceed, CeedVectorNorm(difference_cvec, CEED_NORM_MAX, &max_difference));
45     //TODO Need to reduce across ranks to ensure all ranks are consistent (not a race condition though since the data is purely local)
46     PetscCallCeed(ceed, CeedVectorDestroy(&difference_cvec));
47     if (max_difference < 100 * CEED_EPSILON) {
48       q_data_stored_index = i;
49       break;
50     }
51   }
52 
53   if (q_data_stored_index == -1) {
54     q_data_stored_index = num_q_data_stored++;
55 
56     PetscCall(PetscRealloc(num_q_data_stored * sizeof(CeedVector), &q_data_vecs));
57     PetscCall(PetscRealloc(num_q_data_stored * sizeof(CeedElemRestriction), &q_data_restrictions));
58     q_data_vecs[q_data_stored_index]         = NULL;  // Must set to NULL for ReferenceCopy
59     q_data_restrictions[q_data_stored_index] = NULL;  // Must set to NULL for ReferenceCopy
60     PetscCallCeed(ceed, CeedVectorReferenceCopy(q_data_created, &q_data_vecs[q_data_stored_index]));
61     PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(elem_restr_qd_created, &q_data_restrictions[q_data_stored_index]));
62   }
63   *q_data_stored        = NULL;  // Must set to NULL for ReferenceCopy
64   *elem_restr_qd_stored = NULL;  // Must set to NULL for ReferenceCopy
65   PetscCallCeed(ceed, CeedVectorReferenceCopy(q_data_vecs[q_data_stored_index], q_data_stored));
66   PetscCallCeed(ceed, CeedElemRestrictionReferenceCopy(q_data_restrictions[q_data_stored_index], elem_restr_qd_stored));
67   PetscFunctionReturn(PETSC_SUCCESS);
68 }
69 
70 /**
71   @brief Clear stored QData objects
72 **/
73 PetscErrorCode QDataClearStoredData() {
74   PetscFunctionBeginUser;
75   for (PetscInt i = 0; i < num_q_data_stored; i++) {
76     Ceed ceed = CeedVectorReturnCeed(q_data_vecs[i]);
77 
78     PetscCallCeed(ceed, CeedVectorDestroy(&q_data_vecs[i]));
79     PetscCallCeed(ceed, CeedElemRestrictionDestroy(&q_data_restrictions[i]));
80   }
81   PetscCall(PetscFree(q_data_vecs));
82   PetscCall(PetscFree(q_data_restrictions));
83   q_data_vecs         = NULL;
84   q_data_restrictions = NULL;
85   PetscFunctionReturn(PETSC_SUCCESS);
86 }
87 
88 /**
89  * @brief Get number of components of quadrature data for domain
90  *
91  * @param[in]  dm          DM where quadrature data would be used
92  * @param[out] q_data_size Number of components of quadrature data
93  */
94 PetscErrorCode QDataGetNumComponents(DM dm, CeedInt *q_data_size) {
95   PetscInt num_comp_x, dim;
96 
97   PetscFunctionBeginUser;
98   PetscCall(DMGetDimension(dm, &dim));
99   PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
100   *q_data_size = 0;
101   switch (dim) {
102     case 2:
103       switch (num_comp_x) {
104         case 2:
105           *q_data_size = 5;
106           break;
107         case 3:
108           *q_data_size = 7;
109           break;
110       }
111       break;
112     case 3:
113       *q_data_size = 10;
114       break;
115   }
116   PetscCheck(*q_data_size, 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   PetscFunctionReturn(PETSC_SUCCESS);
119 }
120 
121 /**
122  * @brief Create quadrature data for domain
123  *
124  * @param[in]  ceed          Ceed object quadrature data will be used with
125  * @param[in]  dm            DM where quadrature data would be used
126  * @param[in]  domain_label  DMLabel that quadrature data would be used one
127  * @param[in]  label_value   Value of label
128  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
129  * @param[out] q_data        CeedVector of the quadrature data
130  * @param[out] q_data_size   Number of components of quadrature data
131  */
132 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction *elem_restr_qd, CeedVector *q_data,
133                         CeedInt *q_data_size) {
134   CeedQFunction       qf_setup = NULL;
135   CeedOperator        op_setup;
136   CeedVector          q_data_created, x_coord;
137   CeedElemRestriction elem_restr_qd_created, elem_restr_x;
138   PetscInt            dim, height = 0, num_comp_x;
139   CeedBasis           basis_x;
140 
141   PetscFunctionBeginUser;
142   PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height, &elem_restr_x, &basis_x, &x_coord));
143   PetscCall(QDataGetNumComponents(dm, q_data_size));
144   PetscCall(DMGetCoordinateNumComps(dm, &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, CeedElemRestrictionDestroy(&elem_restr_x));
184   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
185   PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
186   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
187   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
188   PetscFunctionReturn(PETSC_SUCCESS);
189 }
190 
191 /**
192  * @brief Get number of components of quadrature data for boundary of domain
193  *
194  * @param[in]  dm          DM where quadrature data would be used
195  * @param[out] q_data_size Number of components of quadrature data
196  */
197 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) {
198   PetscInt dim;
199 
200   PetscFunctionBeginUser;
201   PetscCall(DMGetDimension(dm, &dim));
202   *q_data_size = 0;
203   switch (dim) {
204     case 2:
205       *q_data_size = 3;
206       break;
207     case 3:
208       *q_data_size = 10;
209       break;
210   }
211   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
212   PetscFunctionReturn(PETSC_SUCCESS);
213 }
214 
215 /**
216  * @brief Create quadrature data for boundary of domain
217  *
218  * @param[in]  ceed          Ceed object quadrature data will be used with
219  * @param[in]  dm            DM where quadrature data would be used
220  * @param[in]  domain_label  DMLabel that quadrature data would be used one
221  * @param[in]  label_value   Value of label
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_qd, CeedVector *q_data,
227                                 CeedInt *q_data_size) {
228   CeedQFunction       qf_setup_sur = NULL;
229   CeedOperator        op_setup_sur;
230   CeedVector          q_data_created, x_coord;
231   CeedElemRestriction elem_restr_qd_created, elem_restr_x;
232   CeedBasis           basis_x;
233   PetscInt            dim, height = 1, num_comp_x;
234 
235   PetscFunctionBeginUser;
236   PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height, &elem_restr_x, &basis_x, &x_coord));
237   PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
238   PetscCall(DMGetCoordinateNumComps(dm, &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, CeedElemRestrictionDestroy(&elem_restr_x));
271   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
272   PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
273   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
274   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
275   PetscFunctionReturn(PETSC_SUCCESS);
276 }
277 
278 /**
279  * @brief Get number of components of quadrature data for boundary of domain
280  *
281  * @param[in]  dm          DM where quadrature data would be used
282  * @param[out] q_data_size Number of components of quadrature data
283  */
284 PetscErrorCode QDataBoundaryGradientGetNumComponents(DM dm, CeedInt *q_data_size) {
285   PetscInt dim;
286 
287   PetscFunctionBeginUser;
288   PetscCall(DMGetDimension(dm, &dim));
289   *q_data_size = 0;
290   switch (dim) {
291     case 2:
292       *q_data_size = 7;
293       break;
294     case 3:
295       *q_data_size = 13;
296       break;
297   }
298   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
299   PetscFunctionReturn(PETSC_SUCCESS);
300 }
301 
302 /**
303   @brief Compute `CeedOperator` surface gradient QData
304 
305   Collective across MPI processes.
306 
307   @param[in]  ceed          `Ceed` object
308   @param[in]  dm            `DMPlex` grid
309   @param[in]  domain_label  `DMLabel` for surface
310   @param[in]  label_value   `DMPlex` label value for surface
311   @param[in]  x_coord       `CeedVector` for coordinates
312   @param[out] elem_restr_qd `CeedElemRestriction` for QData
313   @param[out] q_data        `CeedVector` holding QData
314   @param[out] q_data_size   Number of QData components per quadrature point
315 **/
316 PetscErrorCode QDataBoundaryGradientGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedVector x_coord,
317                                         CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
318   PetscInt            dim;
319   const PetscInt      height_cell = 0, height_face = 1;
320   CeedInt             num_comp_x;
321   CeedElemRestriction elem_restr_x_cell, elem_restr_x_face, elem_restr_qd_created;
322   CeedVector          q_data_created;
323   CeedBasis           basis_x_cell_to_face, basis_x_face;
324   CeedQFunction       qf_setup_sur = NULL;
325   CeedOperator        op_setup_sur;
326 
327   PetscFunctionBeginUser;
328   PetscCall(DMPlexCeedBasisCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &basis_x_face));
329   PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, domain_label, label_value, label_value, &basis_x_cell_to_face));
330   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_cell, &elem_restr_x_cell));
331   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &elem_restr_x_face));
332 
333   PetscCall(QDataBoundaryGradientGetNumComponents(dm, q_data_size));
334   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height_face, *q_data_size, &elem_restr_qd_created));
335   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
336   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x_face, &num_comp_x));
337 
338   PetscCall(DMGetDimension(dm, &dim));
339   switch (dim) {
340     case 2:
341       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2DBoundaryGradient, Setup2DBoundaryGradient_loc, &qf_setup_sur));
342       break;
343     case 3:
344       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundaryGradient, SetupBoundaryGradient_loc, &qf_setup_sur));
345       break;
346   }
347   PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
348 
349   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX cell", num_comp_x * (dim - height_cell), CEED_EVAL_GRAD));
350   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX face", num_comp_x * (dim - height_face), CEED_EVAL_GRAD));
351   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
352   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "q data", *q_data_size, CEED_EVAL_NONE));
353 
354   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_sur));
355   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX cell", elem_restr_x_cell, basis_x_cell_to_face, CEED_VECTOR_ACTIVE));
356   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX face", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE));
357   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE));
358   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "q data", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
359 
360   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
361 
362   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
363 
364   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
365   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
366   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
367   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
368   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
369   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
370   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
371   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
372   PetscFunctionReturn(PETSC_SUCCESS);
373 }
374