1 // Copyright (c) 2017-2026, 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, dim_boundary = 2, dXdx_size = num_comp_x * dim_boundary; 18 CeedVector multiplicity, x_stored, scale_stored, stg_data, dXdx; 19 CeedBasis basis_x_to_q_sur; 20 CeedElemRestriction elem_restr_x_sur, elem_restr_q_sur, 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 PetscInt dm_field = 0, height = 1; 25 26 PetscFunctionBeginUser; 27 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label)); 28 29 { // Basis 30 CeedBasis basis_x_sur, basis_q_sur; 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, dm_field, &basis_q_sur)); 37 PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, label_value, height, dm_field, &basis_x_sur)); 38 39 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_sur, basis_q_sur, &basis_x_to_q_sur)); 40 41 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_sur)); 42 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_sur)); 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", num_comp_x * dim_boundary, 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 Setup QFunction 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 // -- Restrictions 61 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, bc->inflows[i], height, dm_field, &elem_restr_q_sur)); 62 PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, bc->inflows[i], height, &elem_restr_x_sur)); 63 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q_sur, &multiplicity, NULL)); 64 PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr_q_sur, multiplicity)); 65 66 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height, num_comp_x, &elem_restr_x_stored)); 67 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height, 1, &elem_restr_scale)); 68 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height, stg_data_size, &elem_restr_stgdata)); 69 PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, bc->inflows[i], height, dXdx_size, &elem_restr_dXdx)); 70 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL)); 71 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_scale, &scale_stored, NULL)); 72 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_stgdata, &stg_data, NULL)); 73 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_dXdx, &dXdx, NULL)); 74 75 // -- Setup Operator 76 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup)); 77 PetscCallCeed(ceed, CeedOperatorSetName(op_setup, "Precomputed data for strong boundary conditions")); 78 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x", elem_restr_x_sur, basis_x_to_q_sur, CEED_VECTOR_ACTIVE)); 79 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dxdX", elem_restr_x_sur, basis_x_to_q_sur, CEED_VECTOR_ACTIVE)); 80 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "multiplicity", elem_restr_q_sur, CEED_BASIS_NONE, multiplicity)); 81 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x stored", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 82 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored)); 83 PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 84 85 // -- Compute geometric factors 86 PetscCallCeed(ceed, CeedOperatorApply(op_setup, ceed_data->x_coord, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE)); 87 88 // -- Compute STGData 89 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stgdata, NULL, NULL, &op_stgdata)); 90 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 91 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 92 PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 93 94 PetscCallCeed(ceed, CeedOperatorApply(op_stgdata, CEED_VECTOR_NONE, stg_data, CEED_REQUEST_IMMEDIATE)); 95 96 // -- Setup BC Sub Operator 97 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_strongbc, NULL, NULL, &op_strong_bc_sub)); 98 PetscCallCeed(ceed, CeedOperatorSetName(op_strong_bc_sub, "Strong STG")); 99 100 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx)); 101 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored)); 102 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored)); 103 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, stg_data)); 104 PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "q", elem_restr_q_sur, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 105 106 // -- Add to composite operator 107 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_strong_bc, op_strong_bc_sub)); 108 109 PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity)); 110 PetscCallCeed(ceed, CeedVectorDestroy(&x_stored)); 111 PetscCallCeed(ceed, CeedVectorDestroy(&scale_stored)); 112 PetscCallCeed(ceed, CeedVectorDestroy(&stg_data)); 113 PetscCallCeed(ceed, CeedVectorDestroy(&dXdx)); 114 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_sur)); 115 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q_sur)); 116 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_stored)); 117 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_scale)); 118 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_stgdata)); 119 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dXdx)); 120 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc_sub)); 121 PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup)); 122 PetscCallCeed(ceed, CeedOperatorDestroy(&op_stgdata)); 123 } 124 125 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_strong_bc, "solution time", &phys->stg_solution_time_label)); 126 127 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_sur)); 128 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_strongbc)); 129 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stgdata)); 130 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup)); 131 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133 134 PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM, 135 Vec cell_geom_FVM, Vec grad_FVM) { 136 Vec boundary_mask; 137 User user; 138 139 PetscFunctionBeginUser; 140 PetscCall(DMGetApplicationContext(dm, &user)); 141 142 if (user->phys->stg_solution_time_label) { 143 PetscCallCeed(user->ceed, CeedOperatorSetContextDouble(user->op_strong_bc_ctx->op, user->phys->stg_solution_time_label, &time)); 144 } 145 146 // Mask Strong BC entries 147 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 148 PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask)); 149 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 150 151 PetscCall(ApplyAddCeedOperatorLocalToLocal(NULL, Q_loc, user->op_strong_bc_ctx)); 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154 155 PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData problem, SimpleBC bc) { 156 CeedOperator op_strong_bc; 157 158 PetscFunctionBeginUser; 159 { 160 Vec boundary_mask, global_vec; 161 162 PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask)); 163 PetscCall(DMGetGlobalVector(dm, &global_vec)); 164 PetscCall(VecZeroEntries(boundary_mask)); 165 PetscCall(VecSet(global_vec, 1.0)); 166 PetscCall(DMGlobalToLocal(dm, global_vec, INSERT_VALUES, boundary_mask)); 167 PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask)); 168 PetscCall(DMRestoreGlobalVector(dm, &global_vec)); 169 } 170 171 PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_strong_bc)); 172 { 173 PetscBool use_strongstg = PETSC_FALSE; 174 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); 175 176 if (use_strongstg) { 177 PetscCall(SetupStrongSTG_Ceed(ceed, ceed_data, dm, problem, bc, user->phys, op_strong_bc)); 178 } 179 } 180 181 PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_strong_bc, CEED_VECTOR_NONE, NULL, NULL, NULL, &user->op_strong_bc_ctx)); 182 183 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_StrongBCCeed)); 184 PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc)); 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187