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
SetupStrongSTG_Ceed(Ceed ceed,CeedData ceed_data,DM dm,ProblemData problem,SimpleBC bc,Physics phys,CeedOperator op_strong_bc)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
DMPlexInsertBoundaryValues_StrongBCCeed(DM dm,PetscBool insert_essential,Vec Q_loc,PetscReal time,Vec face_geom_FVM,Vec cell_geom_FVM,Vec grad_FVM)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
SetupStrongBC_Ceed(Ceed ceed,CeedData ceed_data,DM dm,User user,ProblemData problem,SimpleBC bc)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