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 { // Get number of coordinate components 98 DM dm_coord; 99 PetscSection section_coord; 100 PetscInt field = 0; // Default field has the coordinates 101 102 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 103 PetscCall(DMGetLocalSection(dm_coord, §ion_coord)); 104 PetscCall(PetscSectionGetFieldComponents(section_coord, field, &num_comp_x)); 105 } 106 switch (dim) { 107 case 2: 108 switch (num_comp_x) { 109 case 2: 110 *q_data_size = 5; 111 break; 112 case 3: 113 *q_data_size = 7; 114 break; 115 default: 116 SETERRQ(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 break; 119 } 120 break; 121 case 3: 122 *q_data_size = 10; 123 break; 124 default: 125 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, 126 "QData not valid for DM of dimension %" PetscInt_FMT " and coordinates with dimension %" PetscInt_FMT, dim, num_comp_x); 127 break; 128 } 129 PetscFunctionReturn(PETSC_SUCCESS); 130 } 131 132 /** 133 * @brief Create quadrature data for domain 134 * 135 * @param[in] ceed Ceed object quadrature data will be used with 136 * @param[in] dm DM where quadrature data would be used 137 * @param[in] domain_label DMLabel that quadrature data would be used one 138 * @param[in] label_value Value of label 139 * @param[in] elem_restr_x CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections) 140 * @param[in] basis_x CeedBasis of the coordinates 141 * @param[in] x_coord CeedVector of the coordinates 142 * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data 143 * @param[out] q_data CeedVector of the quadrature data 144 * @param[out] q_data_size Number of components of quadrature data 145 */ 146 PetscErrorCode QDataGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x, 147 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) { 148 CeedQFunction qf_setup = NULL; 149 CeedOperator op_setup; 150 CeedVector q_data_created; 151 CeedElemRestriction elem_restr_qd_created; 152 CeedInt num_comp_x; 153 PetscInt dim, height = 0; 154 155 PetscFunctionBeginUser; 156 PetscCall(QDataGetNumComponents(dm, q_data_size)); 157 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x)); 158 PetscCall(DMGetDimension(dm, &dim)); 159 switch (dim) { 160 case 2: 161 switch (num_comp_x) { 162 case 2: 163 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2d, Setup2d_loc, &qf_setup)); 164 break; 165 case 3: 166 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup2D_3Dcoords, Setup2D_3Dcoords_loc, &qf_setup)); 167 break; 168 } 169 break; 170 case 3: 171 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Setup, Setup_loc, &qf_setup)); 172 break; 173 } 174 175 // -- Create QFunction for quadrature data 176 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup, 0)); 177 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD)); 178 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT)); 179 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "surface qdata", *q_data_size, CEED_EVAL_NONE)); 180 181 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created)); 182 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL)); 183 184 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup)); 185 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE)); 186 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE)); 187 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 188 189 PetscCallCeed(ceed, CeedOperatorApply(op_setup, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE)); 190 191 PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd)); 192 193 PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created)); 194 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created)); 195 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup)); 196 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup)); 197 PetscFunctionReturn(PETSC_SUCCESS); 198 } 199 200 /** 201 * @brief Get number of components of quadrature data for boundary of domain 202 * 203 * @param[in] dm DM where quadrature data would be used 204 * @param[out] q_data_size Number of components of quadrature data 205 */ 206 PetscErrorCode QDataBoundaryGetNumComponents(DM dm, CeedInt *q_data_size) { 207 PetscInt dim; 208 209 PetscFunctionBeginUser; 210 PetscCall(DMGetDimension(dm, &dim)); 211 switch (dim) { 212 case 2: 213 *q_data_size = 3; 214 break; 215 case 3: 216 *q_data_size = 10; 217 break; 218 default: 219 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "QDataBoundary not valid for DM of dimension %" PetscInt_FMT, dim); 220 break; 221 } 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 /** 226 * @brief Create quadrature data for boundary of domain 227 * 228 * @param[in] ceed Ceed object quadrature data will be used with 229 * @param[in] dm DM where quadrature data would be used 230 * @param[in] domain_label DMLabel that quadrature data would be used one 231 * @param[in] label_value Value of label 232 * @param[in] elem_restr_x CeedElemRestriction of the coordinates (must match `domain_label` and `label_value` selections) 233 * @param[in] basis_x CeedBasis of the coordinates 234 * @param[in] x_coord CeedVector of the coordinates 235 * @param[out] elem_restr_qd CeedElemRestriction of the quadrature data 236 * @param[out] q_data CeedVector of the quadrature data 237 * @param[out] q_data_size Number of components of quadrature data 238 */ 239 PetscErrorCode QDataBoundaryGet(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, CeedElemRestriction elem_restr_x, CeedBasis basis_x, 240 CeedVector x_coord, CeedElemRestriction *elem_restr_qd, CeedVector *q_data, CeedInt *q_data_size) { 241 CeedQFunction qf_setup_sur = NULL; 242 CeedOperator op_setup_sur; 243 CeedVector q_data_created; 244 CeedElemRestriction elem_restr_qd_created; 245 CeedInt num_comp_x; 246 PetscInt dim, height = 1; 247 248 PetscFunctionBeginUser; 249 PetscCall(QDataBoundaryGetNumComponents(dm, q_data_size)); 250 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_x, &num_comp_x)); 251 PetscCall(DMGetDimension(dm, &dim)); 252 switch (dim) { 253 case 2: 254 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary2d, SetupBoundary2d_loc, &qf_setup_sur)); 255 break; 256 case 3: 257 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupBoundary, SetupBoundary_loc, &qf_setup_sur)); 258 break; 259 } 260 261 // -- Create QFunction for quadrature data 262 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_setup_sur, 0)); 263 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "dx", num_comp_x * (dim - height), CEED_EVAL_GRAD)); 264 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup_sur, "weight", 1, CEED_EVAL_WEIGHT)); 265 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup_sur, "surface qdata", *q_data_size, CEED_EVAL_NONE)); 266 267 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, domain_label, label_value, height, *q_data_size, &elem_restr_qd_created)); 268 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_qd_created, &q_data_created, NULL)); 269 270 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup_sur, NULL, NULL, &op_setup_sur)); 271 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE)); 272 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE)); 273 PetscCallCeed(ceed, CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_created, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 274 275 PetscCallCeed(ceed, CeedOperatorApply(op_setup_sur, x_coord, q_data_created, CEED_REQUEST_IMMEDIATE)); 276 277 PetscCall(QDataGetStored(q_data_created, elem_restr_qd_created, q_data, elem_restr_qd)); 278 279 PetscCallCeed(ceed, CeedVectorDestroy(&q_data_created)); 280 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd_created)); 281 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup_sur)); 282 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup_sur)); 283 PetscFunctionReturn(PETSC_SUCCESS); 284 } 285