xref: /honee/src/strong_boundary_conditions.c (revision 0aab724992399e7bc3a3c0143aa36f671b31063e) 
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, 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   MPI_Comm            comm = honee->comm;
23 
24   PetscFunctionBeginUser;
25   PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
26   PetscCall(DMGetDimension(dm, &dim));
27   num_comp_x = dim;
28   dXdx_size  = num_comp_x * (dim - height_cell);
29 
30   {  // Basis
31     CeedBasis basis_x_face, basis_q_face;
32     DM        dm_coord;
33 
34     PetscCall(DMGetCoordinateDM(dm, &dm_coord));
35     DMLabel  label       = NULL;
36     PetscInt label_value = 0;
37     PetscCall(CreateBasisFromPlex(ceed, dm, label, label_value, height_face, dm_field, &basis_q_face));
38     PetscCall(CreateBasisFromPlex(ceed, dm_coord, label, label_value, height_face, dm_field, &basis_x_face));
39 
40     PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_face, basis_q_face, &basis_x_to_q_face));
41 
42     PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face));
43     PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_face));
44   }
45 
46   // Setup QFunction
47   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, SetupStrongBC, SetupStrongBC_loc, &qf_setup));
48   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "x", num_comp_x, CEED_EVAL_INTERP));
49   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "dxdX", dXdx_size, CEED_EVAL_GRAD));
50   PetscCallCeed(ceed, CeedQFunctionAddInput(qf_setup, "multiplicity", num_comp_q, CEED_EVAL_NONE));
51   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "x stored", num_comp_x, CEED_EVAL_NONE));
52   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "scale", 1, CEED_EVAL_NONE));
53   PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_setup, "dXdx", dXdx_size, CEED_EVAL_NONE));
54 
55   // Setup STG QFunctions
56   PetscCall(SetupStrongStg_PreProcessing(ceed, problem, num_comp_x, stg_data_size, dXdx_size, &qf_stgdata));
57   PetscCall(SetupStrongStg_QF(ceed, problem, num_comp_x, num_comp_q, stg_data_size, dXdx_size, &qf_strongbc));
58 
59   // Get face label values for the inflow BC
60   PetscInt        num_inflow_faces;
61   const PetscInt *inflow_faces;
62   {
63     BCDefinition bc_def = NULL;
64     for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
65       BCDefinition bc_def_ = problem->bc_defs[b];
66       const char  *name;
67 
68       PetscCall(BCDefinitionGetInfo(bc_def_, &name, NULL, NULL));
69       if (!strcmp(name, "inflow")) {
70         bc_def = bc_def_;
71         break;
72       }
73     }
74     PetscCheck(bc_def, comm, PETSC_ERR_SUP, "Could not find BCDefintion named 'inflow' for STG setup.");
75     PetscCall(BCDefinitionGetInfo(bc_def, NULL, &num_inflow_faces, &inflow_faces));
76   }
77 
78   // Compute contribution on each boundary face
79   for (CeedInt i = 0; i < num_inflow_faces; i++) {
80     DMLabel  face_orientation_label;
81     PetscInt num_orientations_values, *orientation_values;
82 
83     {
84       char *face_orientation_label_name;
85 
86       PetscCall(DMPlexCreateFaceLabel(dm, inflow_faces[i], &face_orientation_label_name));
87       PetscCall(DMGetLabel(dm, face_orientation_label_name, &face_orientation_label));
88       PetscCall(PetscFree(face_orientation_label_name));
89     }
90     PetscCall(DMLabelCreateGlobalValueArray(dm, face_orientation_label, &num_orientations_values, &orientation_values));
91     for (PetscInt o = 0; o < num_orientations_values; o++) {
92       CeedBasis           basis_x_to_q_cell;
93       CeedElemRestriction elem_restr_x_cell;
94       PetscInt            orientation = orientation_values[o];
95 
96       {
97         CeedBasis basis_x_cell_to_face, basis_q_face;
98 
99         PetscCall(DMPlexCeedBasisCellToFaceCoordinateCreate(ceed, dm, face_orientation_label, orientation, orientation, &basis_x_cell_to_face));
100         PetscCall(CreateBasisFromPlex(ceed, dm, face_orientation_label, orientation, height_face, dm_field, &basis_q_face));
101         PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x_cell_to_face, basis_q_face, &basis_x_to_q_cell));
102         PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_cell_to_face));
103         PetscCallCeed(ceed, CeedBasisDestroy(&basis_q_face));
104       }
105 
106       PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, domain_label, inflow_faces[i], height_face, dm_field, &elem_restr_q_face));
107       PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, domain_label, inflow_faces[i], height_face, &elem_restr_x_face));
108       PetscCall(DMPlexCeedElemRestrictionCoordinateCreate(ceed, dm, face_orientation_label, orientation, height_cell, &elem_restr_x_cell));
109       PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q_face, &multiplicity, NULL));
110       PetscCallCeed(ceed, CeedElemRestrictionGetMultiplicity(elem_restr_q_face, multiplicity));
111 
112       PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, inflow_faces[i], height_face, num_comp_x, &elem_restr_x_stored));
113       PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, inflow_faces[i], height_face, 1, &elem_restr_scale));
114       PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, inflow_faces[i], height_face, stg_data_size, &elem_restr_stgdata));
115       PetscCall(DMPlexCeedElemRestrictionCollocatedCreate(ceed, dm, domain_label, inflow_faces[i], height_face, dXdx_size, &elem_restr_dXdx));
116       PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_x_stored, &x_stored, NULL));
117       PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_scale, &scale_stored, NULL));
118       PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_stgdata, &stg_data, NULL));
119       PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_dXdx, &dXdx, NULL));
120 
121       // -- Setup Operator
122       PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup));
123       PetscCallCeed(ceed, CeedOperatorSetName(op_setup, "Precomputed data for strong boundary conditions"));
124       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x", elem_restr_x_face, basis_x_to_q_face, CEED_VECTOR_ACTIVE));
125       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dxdX", elem_restr_x_cell, basis_x_to_q_cell, CEED_VECTOR_ACTIVE));
126       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "multiplicity", elem_restr_q_face, CEED_BASIS_NONE, multiplicity));
127       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "x stored", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
128       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored));
129       PetscCallCeed(ceed, CeedOperatorSetField(op_setup, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
130 
131       // -- Compute geometric factors
132       PetscCallCeed(ceed, CeedOperatorApply(op_setup, honee->x_coord, CEED_VECTOR_NONE, CEED_REQUEST_IMMEDIATE));
133 
134       // -- Compute STGData
135       PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_stgdata, NULL, NULL, &op_stgdata));
136       PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
137       PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
138       PetscCallCeed(ceed, CeedOperatorSetField(op_stgdata, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
139 
140       PetscCallCeed(ceed, CeedOperatorApply(op_stgdata, CEED_VECTOR_NONE, stg_data, CEED_REQUEST_IMMEDIATE));
141 
142       // -- Setup BC QFunctions
143       PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_strongbc, NULL, NULL, &op_strong_bc_sub));
144       PetscCallCeed(ceed, CeedOperatorSetName(op_strong_bc_sub, "Strong STG"));
145 
146       PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "dXdx", elem_restr_dXdx, CEED_BASIS_NONE, dXdx));
147       PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "x", elem_restr_x_stored, CEED_BASIS_NONE, x_stored));
148       PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "scale", elem_restr_scale, CEED_BASIS_NONE, scale_stored));
149       PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "stg data", elem_restr_stgdata, CEED_BASIS_NONE, stg_data));
150       PetscCallCeed(ceed, CeedOperatorSetField(op_strong_bc_sub, "q", elem_restr_q_face, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
151 
152       // -- Add to composite operator
153       PetscCallCeed(ceed, CeedCompositeOperatorAddSub(op_strong_bc, op_strong_bc_sub));
154 
155       PetscCallCeed(ceed, CeedVectorDestroy(&multiplicity));
156       PetscCallCeed(ceed, CeedVectorDestroy(&x_stored));
157       PetscCallCeed(ceed, CeedVectorDestroy(&scale_stored));
158       PetscCallCeed(ceed, CeedVectorDestroy(&stg_data));
159       PetscCallCeed(ceed, CeedVectorDestroy(&dXdx));
160       PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_cell));
161       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_face));
162       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_cell));
163       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q_face));
164       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x_stored));
165       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_scale));
166       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_stgdata));
167       PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_dXdx));
168       PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc_sub));
169       PetscCallCeed(ceed, CeedOperatorDestroy(&op_setup));
170       PetscCallCeed(ceed, CeedOperatorDestroy(&op_stgdata));
171     }
172     PetscCall(PetscFree(orientation_values));
173   }
174 
175   PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_strong_bc, "solution time", &phys->stg_solution_time_label));
176 
177   PetscCallCeed(ceed, CeedBasisDestroy(&basis_x_to_q_face));
178   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_strongbc));
179   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_stgdata));
180   PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_setup));
181   PetscFunctionReturn(PETSC_SUCCESS);
182 }
183 
184 PetscErrorCode DMPlexInsertBoundaryValues_StrongBCCeed(DM dm, PetscBool insert_essential, Vec Q_loc, PetscReal time, Vec face_geom_FVM,
185                                                        Vec cell_geom_FVM, Vec grad_FVM) {
186   Vec   boundary_mask;
187   Honee honee;
188 
189   PetscFunctionBeginUser;
190   PetscCall(PetscLogEventBegin(HONEE_StrongBCInsert, dm, Q_loc, 0, 0));
191   PetscCall(DMGetApplicationContext(dm, &honee));
192 
193   if (honee->phys->stg_solution_time_label) {
194     PetscCallCeed(honee->ceed, CeedOperatorSetContextDouble(honee->op_strong_bc_ctx->op, honee->phys->stg_solution_time_label, &time));
195   }
196 
197   // Mask Strong BC entries
198   PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask));
199   PetscCall(VecPointwiseMult(Q_loc, Q_loc, boundary_mask));
200   PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask));
201 
202   PetscCall(PetscLogEventBegin(HONEE_StrongBCCeed, dm, Q_loc, 0, 0));
203   PetscCall(ApplyAddCeedOperatorLocalToLocal(NULL, Q_loc, honee->op_strong_bc_ctx));
204   PetscCall(PetscLogEventEnd(HONEE_StrongBCCeed, dm, Q_loc, 0, 0));
205   PetscCall(PetscLogEventEnd(HONEE_StrongBCInsert, dm, Q_loc, 0, 0));
206   PetscFunctionReturn(PETSC_SUCCESS);
207 }
208 
209 PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, DM dm, Honee honee, ProblemData problem) {
210   CeedOperator op_strong_bc;
211 
212   PetscFunctionBeginUser;
213   {
214     Vec boundary_mask, global_vec;
215 
216     PetscCall(DMGetNamedLocalVector(dm, "boundary mask", &boundary_mask));
217     PetscCall(DMGetGlobalVector(dm, &global_vec));
218     PetscCall(VecZeroEntries(boundary_mask));
219     PetscCall(VecSet(global_vec, 1.0));
220     PetscCall(DMGlobalToLocal(dm, global_vec, INSERT_VALUES, boundary_mask));
221     PetscCall(DMRestoreNamedLocalVector(dm, "boundary mask", &boundary_mask));
222     PetscCall(DMRestoreGlobalVector(dm, &global_vec));
223   }
224 
225   PetscCallCeed(ceed, CeedCompositeOperatorCreate(ceed, &op_strong_bc));
226   {
227     PetscBool use_strongstg = PETSC_FALSE;
228     PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL));
229 
230     if (use_strongstg) {
231       PetscCall(SetupStrongSTG_Ceed(ceed, honee, dm, problem, honee->phys, op_strong_bc));
232     }
233   }
234 
235   PetscCall(OperatorApplyContextCreate(NULL, NULL, ceed, op_strong_bc, CEED_VECTOR_NONE, NULL, NULL, NULL, &honee->op_strong_bc_ctx));
236 
237   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertBoundaryValues_C", DMPlexInsertBoundaryValues_StrongBCCeed));
238   PetscCallCeed(ceed, CeedOperatorDestroy(&op_strong_bc));
239   PetscFunctionReturn(PETSC_SUCCESS);
240 }
241