xref: /honee/src/qdata.c (revision 9018c49aeac90a17270a47871327cf384358a98c) 
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[in]  elem_restr_x  CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
223  * @param[in]  basis_x       CeedBasis of the coordinates
224  * @param[in]  x_coord       CeedVector of the coordinates
225  * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
226  * @param[out] q_data        CeedVector of the quadrature data
227  * @param[out] q_data_size   Number of components of quadrature data
228  */
229 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
230                                 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
231   CeedQFunction       qf_setup_sur = NULL;
232   CeedOperator        op_setup_sur;
233   CeedVector          q_data_created;
234   CeedElemRestriction elem_restr_qd_created;
235   CeedInt             num_comp_x;
236   PetscInt            dim, height = 1;
237 
238   PetscFunctionBeginUser;
239   PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
240   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
241   PetscCall(DMGetDimension(dm, &dim));
242   switch (dim) {
243     case 2:
244       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur));
245       break;
246     case 3:
247       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur));
248       break;
249   }
250   PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
251 
252   // -- Create QFunction for quadrature data
253   PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0));
254   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
255   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
256   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE));
257 
258   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created));
259   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
260 
261   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur));
262   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
263   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
264   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
265 
266   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
267 
268   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
269 
270   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
271   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
272   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
273   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
274   PetscFunctionReturn(PETSC_SUCCESS);
275 }
276 
277 /**
278  * @brief Get number of components of quadrature data for boundary of domain
279  *
280  * @param[in]  dm          DM where quadrature data would be used
281  * @param[out] q_data_size Number of components of quadrature data
282  */
283 PetscErrorCode QDataBoundaryGradientGetNumComponents(DM dm, CeedInt *q_data_size) {
284   PetscInt dim;
285 
286   PetscFunctionBeginUser;
287   PetscCall(DMGetDimension(dm, &dim));
288   *q_data_size = 0;
289   switch (dim) {
290     case 2:
291       *q_data_size = 7;
292       break;
293     case 3:
294       *q_data_size = 13;
295       break;
296   }
297   PetscCheck(*q_data_size, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
298   PetscFunctionReturn(PETSC_SUCCESS);
299 }
300 
301 /**
302   @brief Compute `CeedOperator` surface gradient QData
303 
304   Collective across MPI processes.
305 
306   @param[in]  ceed          `Ceed` object
307   @param[in]  dm            `DMPlex` grid
308   @param[in]  domain_label  `DMLabel` for surface
309   @param[in]  label_value   `DMPlex` label value for surface
310   @param[in]  x_coord       `CeedVector` for coordinates
311   @param[out] elem_restr_qd `CeedElemRestriction` for QData
312   @param[out] q_data        `CeedVector` holding QData
313   @param[out] q_data_size   Number of QData components per quadrature point
314 **/
315 PetscErrorCode QDataBoundaryGradientGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedVector x_coord,
316                                         CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
317   PetscInt            dim;
318   const PetscInt      height_cell = 0, height_face = 1;
319   CeedInt             num_comp_x;
320   CeedElemRestriction elem_restr_x_cell, elem_restr_x_face, elem_restr_qd_created;
321   CeedVector          q_data_created;
322   CeedBasis           basis_x_cell_to_face, basis_x_face;
323   CeedQFunction       qf_setup_sur = NULL;
324   CeedOperator        op_setup_sur;
325 
326   PetscFunctionBeginUser;
327   PetscCall(DMPlexCeedBasisCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &basis_x_face));
328   PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, domain_label, label_value, label_value, &basis_x_cell_to_face));
329   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_cell, &elem_restr_x_cell));
330   PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, label_value, height_face, &elem_restr_x_face));
331 
332   PetscCall(QDataBoundaryGradientGetNumComponents(dm, q_data_size));
333   PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height_face, *q_data_size, &elem_restr_qd_created));
334   PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
335   PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x_face, &num_comp_x));
336 
337   PetscCall(DMGetDimension(dm, &dim));
338   switch (dim) {
339     case 2:
340       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2DBoundaryGradient, Setup2DBoundaryGradient_loc, &qf_setup_sur));
341       break;
342     case 3:
343       PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundaryGradient, SetupBoundaryGradient_loc, &qf_setup_sur));
344       break;
345   }
346   PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
347 
348   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX cell", num_comp_x * (dim - height_cell), CEED_EVAL_GRAD));
349   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX face", num_comp_x * (dim - height_face), CEED_EVAL_GRAD));
350   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
351   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "q data", *q_data_size, CEED_EVAL_NONE));
352 
353   PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_sur));
354   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX cell", elem_restr_x_cell, basis_x_cell_to_face, CEED_VECTOR_ACTIVE));
355   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX face", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE));
356   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE));
357   PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "q data", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
358 
359   PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
360 
361   PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
362 
363   PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
364   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
365   PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
366   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created));
367   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
368   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
369   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
370   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
371   PetscFunctionReturn(PETSC_SUCCESS);
372 }
373