1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #include "../navierstokes.h"
9
10 #include <petscsection.h>
11 #include "../qfunctions/setupgeo.h"
12 #include "../qfunctions/setupgeo2d.h"
13
14 /**
15 * @brief Get number of components of quadrature data for domain
16 *
17 * @param[in] dm DM where quadrature data would be used
18 * @param[out] q_data_size Number of components of quadrature data
19 */
QDataGetNumComponents(DM dm,CeedInt * q_data_size)20 PetscErrorCode QDataGetNumComponents(DM dm, CeedInt *q_data_size) {
21 PetscInt num_comp_x, dim;
22
23 PetscFunctionBeginUser;
24 PetscCall(DMGetDimension(dm, &dim));
25 { // Get number of coordinate components
26 DM dm_coord;
27 PetscSection section_coord;
28 PetscInt field = 0; // Default field has the coordinates
29 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
30 PetscCall(DMGetLocalSection(dm_coord, §ion_coord));
31 PetscCall(PetscSectionGetFieldComponents(section_coord, field, &num_comp_x));
32 }
33 switch (dim) {
34 case 2:
35 switch (num_comp_x) {
36 case 2:
37 *q_data_size = 5;
38 break;
39 case 3:
40 *q_data_size = 7;
41 break;
42 default:
43 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
44 "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
45 break;
46 }
47 break;
48 case 3:
49 *q_data_size = 10;
50 break;
51 default:
52 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,
53 "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x);
54 break;
55 }
56 PetscFunctionReturn(PETSC_SUCCESS);
57 }
58
59 /**
60 * @brief Create quadrature data for domain
61 *
62 * @param[in] ceed Ceed object quadrature data will be used with
63 * @param[in] dm DM where quadrature data would be used
64 * @param[in] domain_label DMLabel that quadrature data would be used one
65 * @param[in] label_value Value of label
66 * @param[in] elem_restr_x CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
67 * @param[in] basis_x CeedBasis of the coordinates
68 * @param[in] x_coord CeedVector of the coordinates
69 * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
70 * @param[out] q_data CeedVector of the quadrature data
71 * @param[out] q_data_size number of components of quadrature data
72 */
QDataGet(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,CeedElemRestriction elem_restr_x,CeedBasis basis_x,CeedVector x_coord,CeedElemRestriction * elem_restr_qd,CeedVector * q_data,CeedInt * q_data_size)73 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
74 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
75 CeedQFunction qf_setup;
76 CeedOperator op_setup;
77 CeedInt num_comp_x;
78 PetscInt dim, height = 0;
79
80 PetscFunctionBeginUser;
81 PetscCall(QDataGetNumComponents(dm, q_data_size));
82 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
83 PetscCall(DMGetDimension(dm, &dim));
84 switch (dim) {
85 case 2:
86 switch (num_comp_x) {
87 case 2:
88 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2d, Setup2d_loc, &qf_setup));
89 break;
90 case 3:
91 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2D_3Dcoords, Setup2D_3Dcoords_loc, &qf_setup));
92 break;
93 }
94 break;
95 case 3:
96 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup, Setup_loc, &qf_setup));
97 break;
98 }
99
100 // -- Create QFunction for quadrature data
101 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup, 0));
102 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
103 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT));
104 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "surface qdata", *q_data_size, CEED_EVAL_NONE));
105
106 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, elem_restr_qd));
107 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_qd, q_data, NULL));
108
109 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
110 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
111 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
112 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "surface qdata", *elem_restr_qd, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
113
114 PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, *q_data, CEED_REQUEST_IMMEDIATE));
115
116 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
117 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
118 PetscFunctionReturn(PETSC_SUCCESS);
119 }
120
121 /**
122 * @brief Get number of components of quadrature data for boundary of domain
123 *
124 * @param[in] dm DM where quadrature data would be used
125 * @param[out] q_data_size Number of components of quadrature data
126 */
QDataBoundaryGetNumComponents(DM dm,CeedInt * q_data_size)127 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) {
128 PetscInt dim;
129
130 PetscFunctionBeginUser;
131 PetscCall(DMGetDimension(dm, &dim));
132 switch (dim) {
133 case 2:
134 *q_data_size = 3;
135 break;
136 case 3:
137 *q_data_size = 10;
138 break;
139 default:
140 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "QDataBoundary not valid for DM of dimension %" PetscInt_FMT, dim);
141 break;
142 }
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
146 /**
147 * @brief Create quadrature data for boundary of domain
148 *
149 * @param[in] ceed Ceed object quadrature data will be used with
150 * @param[in] dm DM where quadrature data would be used
151 * @param[in] domain_label DMLabel that quadrature data would be used one
152 * @param[in] label_value Value of label
153 * @param[in] elem_restr_x CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections)
154 * @param[in] basis_x CeedBasis of the coordinates
155 * @param[in] x_coord CeedVector of the coordinates
156 * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data
157 * @param[out] q_data CeedVector of the quadrature data
158 * @param[out] q_data_size number of components of quadrature data
159 */
QDataBoundaryGet(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,CeedElemRestriction elem_restr_x,CeedBasis basis_x,CeedVector x_coord,CeedElemRestriction * elem_restr_qd,CeedVector * q_data,CeedInt * q_data_size)160 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x,
161 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) {
162 CeedQFunction qf_setup_sur;
163 CeedOperator op_setup_sur;
164 CeedInt num_comp_x;
165 PetscInt dim, height = 1;
166
167 PetscFunctionBeginUser;
168 PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size));
169 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x));
170 PetscCall(DMGetDimension(dm, &dim));
171 switch (dim) {
172 case 2:
173 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur));
174 break;
175 case 3:
176 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur));
177 break;
178 }
179
180 // -- Create QFunction for quadrature data
181 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0));
182 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD));
183 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT));
184 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE));
185
186 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, elem_restr_qd));
187 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(*elem_restr_qd, q_data, NULL));
188
189 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur));
190 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE));
191 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE));
192 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", *elem_restr_qd, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
193
194 PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, *q_data, CEED_REQUEST_IMMEDIATE));
195
196 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur));
197 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur));
198 PetscFunctionReturn(PETSC_SUCCESS);
199 }
200