1 // Copyright (c) 2017-2024, 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 */ 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 */ 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 */ 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 */ 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