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 **/
QDataGetStored(CeedVector q_data_created,CeedElemRestriction elem_restr_qd_created,CeedVector * q_data_stored,CeedElemRestriction * elem_restr_qd_stored)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 **/
QDataClearStoredData()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 */
QDataGetNumComponents(DM dm,CeedInt * q_data_size)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 */
QDataGet(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,CeedElemRestriction * elem_restr_qd,CeedVector * q_data,CeedInt * q_data_size)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 */
QDataBoundaryGetNumComponents(DM dm,CeedInt * q_data_size)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 */
QDataBoundaryGet(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,CeedElemRestriction * elem_restr_qd,CeedVector * q_data,CeedInt * q_data_size)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 */
QDataBoundaryGradientGetNumComponents(DM dm,CeedInt * q_data_size)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[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 **/
QDataBoundaryGradientGet(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,CeedElemRestriction * elem_restr_qd,CeedVector * q_data,CeedInt * q_data_size)315 PetscErrorCode QDataBoundaryGradientGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction *elem_restr_qd,
316 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, x_coord;
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(DMPlexCeedCoordinateCreateField(ceed, dm, domain_label, label_value, height_face, &elem_restr_x_face, &basis_x_face, &x_coord));
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
331 PetscCall(QDataBoundaryGradientGetNumComponents(dm, q_data_size));
332 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height_face, *q_data_size, &elem_restr_qd_created));
333 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL));
334 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x_face, &num_comp_x));
335
336 PetscCall(DMGetDimension(dm, &dim));
337 switch (dim) {
338 case 2:
339 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2DBoundaryGradient, Setup2DBoundaryGradient_loc, &qf_setup_sur));
340 break;
341 case 3:
342 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundaryGradient, SetupBoundaryGradient_loc, &qf_setup_sur));
343 break;
344 }
345 PetscCheck(qf_setup_sur, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "%s not valid for DM of dimension %" PetscInt_FMT, __func__, dim);
346
347 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX cell", num_comp_x * (dim - height_cell), CEED_EVAL_GRAD));
348 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx/dX face", num_comp_x * (dim - height_face), CEED_EVAL_GRAD));
349 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
350 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "q data", *q_data_size, CEED_EVAL_NONE));
351
352 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_sur));
353 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX cell", elem_restr_x_cell, basis_x_cell_to_face, CEED_VECTOR_ACTIVE));
354 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx/dX face", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE));
355 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE));
356 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "q data", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
357
358 PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE));
359
360 PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd));
361
362 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
363 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
364 PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created));
365 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
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