1 // Copyright (c) 2017-2022, 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, CeedInt q_data_size_sur, 17 CeedOperator op_strong_bc) { 18 CeedInt num_comp_x = problem->dim, num_comp_q = 5, num_elem, elem_size, stg_data_size = 1; 19 CeedVector multiplicity, x_stored, scale_stored, q_data_sur, stg_data; 20 CeedBasis basis_x_to_q_sur; 21 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, elem_restr_x_stored, elem_restr_scale, elem_restr_qd_sur, elem_restr_stgdata; 22 CeedQFunction qf_setup, qf_strongbc, qf_stgdata; 23 CeedOperator op_setup, op_strong_bc_sub, op_setup_sur, op_stgdata; 24 DMLabel domain_label; 25 26 PetscFunctionBeginUser; 27 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label)); 28 29 // Basis 30 CeedInt height = 1; 31 PetscCall(CeedBasisCreateProjection(ceed_data->basis_x_sur, ceed_data->basis_q_sur, &basis_x_to_q_sur)); 32 // --- Get number of quadrature points for the boundaries 33 CeedInt num_qpts_sur; 34 CeedBasisGetNumQuadraturePoints(ceed_data->basis_q_sur, &num_qpts_sur); 35 36 // Setup QFunction 37 CeedQFunctionCreateInterior(ceed, 1, SetupStrongBC, SetupStrongBC_loc, &qf_setup); 38 CeedQFunctionAddInput(qf_setup, "x", num_comp_x, CEED_EVAL_INTERP); 39 CeedQFunctionAddInput(qf_setup, "multiplicity", num_comp_q, CEED_EVAL_NONE); 40 CeedQFunctionAddOutput(qf_setup, "x stored", num_comp_x, CEED_EVAL_NONE); 41 CeedQFunctionAddOutput(qf_setup, "scale", 1, CEED_EVAL_NONE); 42 43 // Setup STG Setup QFunction 44 PetscCall(SetupStrongSTG_PreProcessing(ceed, problem, num_comp_x, stg_data_size, q_data_size_sur, &qf_stgdata)); 45 46 // Compute contribution on each boundary face 47 for (CeedInt i = 0; i < bc->num_inflow; i++) { 48 // -- Restrictions 49 PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, -1, -1, &elem_restr_q_sur, &elem_restr_x_sur, NULL)); 50 CeedElemRestrictionCreateVector(elem_restr_q_sur, &multiplicity, NULL); 51 CeedElemRestrictionGetMultiplicity(elem_restr_q_sur, multiplicity); 52 CeedElemRestrictionGetNumElements(elem_restr_q_sur, &num_elem); 53 CeedElemRestrictionGetElementSize(elem_restr_q_sur, &elem_size); 54 PetscCall(GetRestrictionForDomain(ceed, dm, height, domain_label, bc->inflows[i], 0, elem_size, q_data_size_sur, NULL, NULL, &elem_restr_qd_sur)); 55 56 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, num_comp_x, num_elem * elem_size * num_comp_x, CEED_STRIDES_BACKEND, 57 &elem_restr_x_stored); 58 CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL); 59 60 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, 1, num_elem * elem_size, CEED_STRIDES_BACKEND, &elem_restr_scale); 61 CeedElemRestrictionCreateVector(elem_restr_scale, &scale_stored, NULL); 62 63 CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, stg_data_size, num_elem * elem_size, CEED_STRIDES_BACKEND, &elem_restr_stgdata); 64 CeedElemRestrictionCreateVector(elem_restr_stgdata, &stg_data, NULL); 65 66 CeedVectorCreate(ceed, q_data_size_sur * num_elem * elem_size, &q_data_sur); 67 68 // -- Setup Operator 69 CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup); 70 CeedOperatorSetName(op_setup, "surface geometric data"); 71 CeedOperatorSetField(op_setup, "x", elem_restr_x_sur, basis_x_to_q_sur, CEED_VECTOR_ACTIVE); 72 CeedOperatorSetField(op_setup, "multiplicity", elem_restr_q_sur, CEED_BASIS_COLLOCATED, multiplicity); 73 CeedOperatorSetField(op_setup, "x stored", elem_restr_x_stored, CEED_BASIS_COLLOCATED, x_stored); 74 CeedOperatorSetField(op_setup, "scale", elem_restr_scale, CEED_BASIS_COLLOCATED, scale_stored); 75 76 // -- Compute geometric factors 77 CeedOperatorApply(op_setup, ceed_data->x_coord, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE); 78 79 // -- Compute QData for the surface 80 CeedOperatorCreate(ceed, ceed_data->qf_setup_sur, NULL, NULL, &op_setup_sur); 81 CeedOperatorSetField(op_setup_sur, "dx", elem_restr_x_sur, ceed_data->basis_xc_sur, CEED_VECTOR_ACTIVE); 82 CeedOperatorSetField(op_setup_sur, "weight", CEED_ELEMRESTRICTION_NONE, ceed_data->basis_xc_sur, CEED_VECTOR_NONE); 83 CeedOperatorSetField(op_setup_sur, "surface qdata", elem_restr_qd_sur, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 84 85 CeedOperatorApply(op_setup_sur, ceed_data->x_coord, q_data_sur, CEED_REQUEST_IMMEDIATE); 86 87 // -- Compute STGData 88 CeedOperatorCreate(ceed, qf_stgdata, NULL, NULL, &op_stgdata); 89 CeedOperatorSetField(op_stgdata, "surface qdata", elem_restr_qd_sur, CEED_BASIS_COLLOCATED, q_data_sur); 90 CeedOperatorSetField(op_stgdata, "x", elem_restr_x_stored, CEED_BASIS_COLLOCATED, x_stored); 91 CeedOperatorSetField(op_stgdata, "stg data", elem_restr_stgdata, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 92 93 CeedOperatorApply(op_stgdata, CEED_VECTOR_NONE, stg_data, CEED_REQUEST_IMMEDIATE); 94 95 // -- Setup BC QFunctions 96 SetupStrongSTG_QF(ceed, problem, num_comp_x, num_comp_q, stg_data_size, q_data_size_sur, &qf_strongbc); 97 CeedOperatorCreate(ceed, qf_strongbc, NULL, NULL, &op_strong_bc_sub); 98 CeedOperatorSetName(op_strong_bc_sub, "Strong STG"); 99 100 CeedOperatorSetField(op_strong_bc_sub, "surface qdata", elem_restr_qd_sur, CEED_BASIS_COLLOCATED, q_data_sur); 101 CeedOperatorSetField(op_strong_bc_sub, "x", elem_restr_x_stored, CEED_BASIS_COLLOCATED, x_stored); 102 CeedOperatorSetField(op_strong_bc_sub, "scale", elem_restr_scale, CEED_BASIS_COLLOCATED, scale_stored); 103 CeedOperatorSetField(op_strong_bc_sub, "stg data", elem_restr_stgdata, CEED_BASIS_COLLOCATED, stg_data); 104 CeedOperatorSetField(op_strong_bc_sub, "q", elem_restr_q_sur, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 105 106 // -- Add to composite operator 107 CeedCompositeOperatorAddSub(op_strong_bc, op_strong_bc_sub); 108 109 CeedVectorDestroy(&q_data_sur); 110 CeedVectorDestroy(&multiplicity); 111 CeedVectorDestroy(&x_stored); 112 CeedVectorDestroy(&scale_stored); 113 CeedVectorDestroy(&stg_data); 114 CeedElemRestrictionDestroy(&elem_restr_x_sur); 115 CeedElemRestrictionDestroy(&elem_restr_q_sur); 116 CeedElemRestrictionDestroy(&elem_restr_qd_sur); 117 CeedElemRestrictionDestroy(&elem_restr_x_stored); 118 CeedElemRestrictionDestroy(&elem_restr_scale); 119 CeedElemRestrictionDestroy(&elem_restr_stgdata); 120 CeedQFunctionDestroy(&qf_strongbc); 121 CeedQFunctionDestroy(&qf_stgdata); 122 CeedOperatorDestroy(&op_setup_sur); 123 CeedOperatorDestroy(&op_strong_bc_sub); 124 CeedOperatorDestroy(&op_setup); 125 CeedOperatorDestroy(&op_stgdata); 126 } 127 128 CeedOperatorGetContextFieldLabel(op_strong_bc, "solution time", &phys->stg_solution_time_label); 129 130 CeedBasisDestroy(&basis_x_to_q_sur); 131 CeedQFunctionDestroy(&qf_setup); 132 133 PetscFunctionReturn(PETSC_SUCCESS); 134 } 135 136 PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 137 Vec cell_geom_FVM, Vec grad_FVM) { 138 Vec boundary_mask; 139 User user; 140 141 PetscFunctionBeginUser; 142 PetscCall(DMGetApplicationContext(dm, &user)); 143 144 if (user->phys->stg_solution_time_label) { 145 CeedOperatorSetContextDouble(user->op_strong_bc_ctx->op, user->phys->stg_solution_time_label, &time); 146 } 147 148 // Mask Strong BC entries 149 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 150 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 151 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 152 153 PetscCall(ApplyAddCeedOperatorLocalToLocal(NULL, Q_loc, user->op_strong_bc_ctx)); 154 155 PetscFunctionReturn(PETSC_SUCCESS); 156 } 157 158 PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt q_data_size_sur) { 159 CeedOperator op_strong_bc; 160 161 PetscFunctionBeginUser; 162 { 163 Vec boundary_mask, global_vec; 164 165 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 166 PetscCall(DMGetGlobalVector(dm, &global_vec)); 167 PetscCall(VecZeroEntries(boundary_mask)); 168 PetscCall(VecSet(global_vec, 1.0)); 169 PetscCall(DMGlobalToLocal(dm, global_vec, INSERT_VALUES, boundary_mask)); 170 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 171 PetscCall(DMRestoreGlobalVector(dm, &global_vec)); 172 } 173 174 CeedCompositeOperatorCreate(ceed, &op_strong_bc); 175 { 176 PetscBool use_strongstg = PETSC_FALSE; 177 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 178 179 if (use_strongstg) { 180 PetscCall(SetupStrongSTG_Ceed(ceed, ceed_data, dm, problem, bc, user->phys, q_data_size_sur, op_strong_bc)); 181 } 182 } 183 184 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_strong_bc, CEED_VECTOR_NONE, NULL, NULL, NULL, &user->op_strong_bc_ctx)); 185 186 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_StrongBCCeed)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189