1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include "../qfunctions/strong_boundary_conditions.h" 5 6 #include <ceed.h> 7 #include <petscdmplex.h> 8 9 #include <navierstokes.h> 10 #include "../problems/stg_shur14.h" 11 12 PetscErrorCode SetupStrongSTG_Ceed(Ceed ceed, Honee honee, DM dm, ProblemData problem, SimpleBC bc, Physics phys, CeedOperator op_strong_bc) { 13 CeedInt num_comp_x, num_comp_q = 5, stg_data_size = 1, dXdx_size; 14 CeedVector multiplicity, x_stored, scale_stored, stg_data, dXdx; 15 CeedBasis basis_x_to_q_face; 16 CeedElemRestriction elem_restr_x_face, elem_restr_q_face, elem_restr_x_stored, elem_restr_scale, elem_restr_stgdata, elem_restr_dXdx; 17 CeedQFunction qf_setup, qf_strongbc, qf_stgdata; 18 CeedOperator op_setup, op_strong_bc_sub, op_stgdata; 19 DMLabel domain_label; 20 const PetscInt dm_field = 0, height_face = 1, height_cell = 0; 21 PetscInt dim; 22 23 PetscFunctionBeginUser; 24 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label)); 25 PetscCall(DMGetDimension(dm, &dim)); 26 num_comp_x = dim; 27 dXdx_size = num_comp_x * (dim - height_cell); 28 29 { // Basis 30 CeedBasis basis_x_face, basis_q_face; 31 DM dm_coord; 32 33 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 34 DMLabel label = NULL; 35 PetscInt label_value = 0; 36 PetscCall(CreateBasisFromPlex(ceed, dm, label, label_value, height_face, dm_field, &basis_q_face)); 37 PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, label_value, height_face, dm_field, &basis_x_face)); 38 39 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_face, basis_q_face, &basis_x_to_q_face)); 40 41 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face)); 42 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face)); 43 } 44 45 // Setup QFunction 46 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupStrongBC, SetupStrongBC_loc, &qf_setup)); 47 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "x", num_comp_x, CEED_EVAL_INTERP)); 48 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dxdX", dXdx_size, CEED_EVAL_GRAD)); 49 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "multiplicity", num_comp_q, CEED_EVAL_NONE)); 50 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "x stored", num_comp_x, CEED_EVAL_NONE)); 51 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "scale", 1, CEED_EVAL_NONE)); 52 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "dXdx", dXdx_size, CEED_EVAL_NONE)); 53 54 // Setup STG QFunctions 55 PetscCall(SetupStrongStg_PreProcessing(ceed, problem, num_comp_x, stg_data_size, dXdx_size, &qf_stgdata)); 56 PetscCall(SetupStrongStg_QF(ceed, problem, num_comp_x, num_comp_q, stg_data_size, dXdx_size, &qf_strongbc)); 57 58 // Compute contribution on each boundary face 59 for (CeedInt i = 0; i < bc->num_inflow; i++) { 60 DMLabel face_orientation_label; 61 PetscInt num_orientations_values, *orientation_values; 62 63 { 64 char *face_orientation_label_name; 65 66 PetscCall(DMPlexCreateFaceLabel(dm, bc->inflows[i], &face_orientation_label_name)); 67 PetscCall(DMGetLabel(dm, face_orientation_label_name, &face_orientation_label)); 68 PetscCall(PetscFree(face_orientation_label_name)); 69 } 70 PetscCall(DMLabelCreateGlobalValueArray(dm, face_orientation_label, &num_orientations_values, &orientation_values)); 71 for (PetscInt o = 0; o < num_orientations_values; o++) { 72 CeedBasis basis_x_to_q_cell; 73 CeedElemRestriction elem_restr_x_cell; 74 PetscInt orientation = orientation_values[o]; 75 76 { 77 CeedBasis basis_x_cell_to_face, basis_q_face; 78 79 PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, face_orientation_label, orientation, orientation, &basis_x_cell_to_face)); 80 PetscCall(CreateBasisFromPlex(ceed, dm, face_orientation_label, orientation, height_face, dm_field, &basis_q_face)); 81 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_cell_to_face, basis_q_face, &basis_x_to_q_cell)); 82 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face)); 83 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face)); 84 } 85 86 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, bc->inflows[i], height_face, dm_field, &elem_restr_q_face)); 87 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, bc->inflows[i], height_face, &elem_restr_x_face)); 88 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, face_orientation_label, orientation, height_cell, &elem_restr_x_cell)); 89 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q_face, &multiplicity, NULL)); 90 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr_q_face, multiplicity)); 91 92 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height_face, num_comp_x, &elem_restr_x_stored)); 93 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height_face, 1, &elem_restr_scale)); 94 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height_face, stg_data_size, &elem_restr_stgdata)); 95 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height_face, dXdx_size, &elem_restr_dXdx)); 96 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL)); 97 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_scale, &scale_stored, NULL)); 98 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_stgdata, &stg_data, NULL)); 99 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_dXdx, &dXdx, NULL)); 100 101 // -- Setup Operator 102 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup)); 103 PetscCallCeed(ceed, CeedOperatorSetName(op_setup, "Precomputed data for strong boundary conditions")); 104 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x", elem_restr_x_face, basis_x_to_q_face, CEED_VECTOR_ACTIVE)); 105 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dxdX", elem_restr_x_cell, basis_x_to_q_cell, CEED_VECTOR_ACTIVE)); 106 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "multiplicity", elem_restr_q_face, CEED_BASIS_NONE, multiplicity)); 107 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x stored", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 108 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored)); 109 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 110 111 // -- Compute geometric factors 112 PetscCallCeed(ceed, CeedOperatorApply(op_setup, honee->x_coord, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE)); 113 114 // -- Compute STGData 115 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stgdata, NULL, NULL, &op_stgdata)); 116 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 117 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 118 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 119 120 PetscCallCeed(ceed, CeedOperatorApply(op_stgdata, CEED_VECTOR_NONE, stg_data, CEED_REQUEST_IMMEDIATE)); 121 122 // -- Setup BC QFunctions 123 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_strongbc, NULL, NULL, &op_strong_bc_sub)); 124 PetscCallCeed(ceed, CeedOperatorSetName(op_strong_bc_sub, "Strong STG")); 125 126 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 127 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 128 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored)); 129 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, stg_data)); 130 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "q", elem_restr_q_face, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 131 132 // -- Add to composite operator 133 PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_strong_bc, op_strong_bc_sub)); 134 135 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity)); 136 PetscCallCeed(ceed, CeedVectorDestroy(&x_stored)); 137 PetscCallCeed(ceed, CeedVectorDestroy(&scale_stored)); 138 PetscCallCeed(ceed, CeedVectorDestroy(&stg_data)); 139 PetscCallCeed(ceed, CeedVectorDestroy(&dXdx)); 140 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_cell)); 141 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face)); 142 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell)); 143 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q_face)); 144 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_stored)); 145 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_scale)); 146 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_stgdata)); 147 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dXdx)); 148 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc_sub)); 149 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup)); 150 PetscCallCeed(ceed, CeedOperatorDestroy(&op_stgdata)); 151 } 152 PetscCall(PetscFree(orientation_values)); 153 } 154 155 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_strong_bc, "solution time", &phys->stg_solution_time_label)); 156 157 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_face)); 158 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_strongbc)); 159 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stgdata)); 160 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup)); 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 165 Vec cell_geom_FVM, Vec grad_FVM) { 166 Vec boundary_mask; 167 Honee honee; 168 169 PetscFunctionBeginUser; 170 PetscCall(PetscLogEventBegin(HONEE_StrongBCInsert, dm, Q_loc, 0, 0)); 171 PetscCall(DMGetApplicationContext(dm, &honee)); 172 173 if (honee->phys->stg_solution_time_label) { 174 PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(honee->op_strong_bc_ctx->op, honee->phys->stg_solution_time_label, &time)); 175 } 176 177 // Mask Strong BC entries 178 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 179 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 180 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 181 182 PetscCall(PetscLogEventBegin(HONEE_StrongBCCeed, dm, Q_loc, 0, 0)); 183 PetscCall(ApplyAddCeedOperatorLocalToLocal(NULL, Q_loc, honee->op_strong_bc_ctx)); 184 PetscCall(PetscLogEventEnd(HONEE_StrongBCCeed, dm, Q_loc, 0, 0)); 185 PetscCall(PetscLogEventEnd(HONEE_StrongBCInsert, dm, Q_loc, 0, 0)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, DM dm, Honee honee, ProblemData problem, SimpleBC bc) { 190 CeedOperator op_strong_bc; 191 192 PetscFunctionBeginUser; 193 { 194 Vec boundary_mask, global_vec; 195 196 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 197 PetscCall(DMGetGlobalVector(dm, &global_vec)); 198 PetscCall(VecZeroEntries(boundary_mask)); 199 PetscCall(VecSet(global_vec, 1.0)); 200 PetscCall(DMGlobalToLocal(dm, global_vec, INSERT_VALUES, boundary_mask)); 201 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 202 PetscCall(DMRestoreGlobalVector(dm, &global_vec)); 203 } 204 205 PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, &op_strong_bc)); 206 { 207 PetscBool use_strongstg = PETSC_FALSE; 208 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 209 210 if (use_strongstg) { 211 PetscCall(SetupStrongSTG_Ceed(ceed, honee, dm, problem, bc, honee->phys, op_strong_bc)); 212 } 213 } 214 215 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_strong_bc, CEED_VECTOR_NONE, NULL, NULL, NULL, &honee->op_strong_bc_ctx)); 216 217 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_StrongBCCeed)); 218 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc)); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221