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