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