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