1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Setup libCEED Operators for HONEE
6
7 #include <ceed.h>
8 #include <petscdmplex.h>
9 #include <smartsim.h>
10
11 #include <navierstokes.h>
12
13 // @brief Create CeedOperator for unstabilized mass KSP for explicit timestepping
CreateKSPMassOperator_Unstabilized(Honee honee,CeedOperator * op_mass)14 static PetscErrorCode CreateKSPMassOperator_Unstabilized(Honee honee, CeedOperator *op_mass) {
15 Ceed ceed = honee->ceed;
16 CeedInt num_comp_q, q_data_size;
17 CeedQFunction qf_mass;
18 CeedElemRestriction elem_restr_q, elem_restr_qd;
19 CeedBasis basis_q;
20 CeedVector q_data;
21
22 PetscFunctionBeginUser;
23 { // Get restriction and basis from the RHS function
24 CeedOperator *sub_ops;
25 CeedOperatorField op_field;
26 PetscInt sub_op_index = 0; // will be 0 for the volume op
27
28 PetscCallCeed(ceed, CeedOperatorCompositeGetSubList(honee->op_rhs_ctx->op, &sub_ops));
29 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "q", &op_field));
30 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_q, &basis_q, NULL));
31 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "qdata", &op_field));
32 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_qd, NULL, &q_data));
33 }
34
35 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_q, &num_comp_q));
36 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_qd, &q_data_size));
37
38 PetscCall(HoneeMassQFunctionCreate(ceed, num_comp_q, q_data_size, &qf_mass));
39 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_mass, NULL, NULL, op_mass));
40 PetscCallCeed(ceed, CeedOperatorSetName(*op_mass, "RHS Mass Operator"));
41 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "u", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
42 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
43 PetscCallCeed(ceed, CeedOperatorSetField(*op_mass, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
44
45 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
46 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
47 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
48 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
49 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_mass));
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
53 // @brief Create KSP to solve the inverse mass operator for explicit time stepping schemes
CreateKSPMass(Honee honee,ProblemData problem)54 static PetscErrorCode CreateKSPMass(Honee honee, ProblemData problem) {
55 Ceed ceed = honee->ceed;
56 DM dm = honee->dm;
57 CeedOperator op_mass;
58
59 PetscFunctionBeginUser;
60 if (problem->create_mass_operator) PetscCall(problem->create_mass_operator(honee, &op_mass));
61 else PetscCall(CreateKSPMassOperator_Unstabilized(honee, &op_mass));
62
63 { // -- Setup KSP for mass operator
64 Mat mat_mass;
65 Vec Zeros_loc;
66 MPI_Comm comm = PetscObjectComm((PetscObject)dm);
67
68 PetscCall(DMCreateLocalVector(dm, &Zeros_loc));
69 PetscCall(VecZeroEntries(Zeros_loc));
70 PetscCall(MatCreateCeed(dm, dm, op_mass, NULL, &mat_mass));
71 PetscCall(MatCeedSetLocalVectors(mat_mass, Zeros_loc, NULL));
72
73 PetscCall(KSPCreate(comm, &honee->mass_ksp));
74 PetscCall(KSPSetOptionsPrefix(honee->mass_ksp, "mass_"));
75 PetscCall(PetscObjectSetName((PetscObject)honee->mass_ksp, "Explicit Mass"));
76 { // lumped by default
77 PC pc;
78 PetscCall(KSPGetPC(honee->mass_ksp, &pc));
79 PetscCall(PCSetType(pc, PCJACOBI));
80 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
81 PetscCall(KSPSetType(honee->mass_ksp, KSPPREONLY));
82 }
83 PetscCall(KSPSetFromOptions_WithMatCeed(honee->mass_ksp, mat_mass));
84 PetscCall(VecDestroy(&Zeros_loc));
85 PetscCall(MatDestroy(&mat_mass));
86 }
87
88 PetscCallCeed(ceed, CeedOperatorDestroy(&op_mass));
89 PetscFunctionReturn(PETSC_SUCCESS);
90 }
91
SetupLibceed(Ceed ceed,DM dm,Honee honee,AppCtx app_ctx,ProblemData problem)92 PetscErrorCode SetupLibceed(Ceed ceed, DM dm, Honee honee, AppCtx app_ctx, ProblemData problem) {
93 const PetscInt num_comp_q = problem->num_components;
94 PetscInt dim, num_comp_x;
95 CeedInt q_data_size_vol, num_comps_jac_data = problem->num_comps_jac_data;
96 CeedElemRestriction elem_restr_qd, elem_restr_q, elem_restr_x, elem_restr_jd_i = NULL;
97 CeedBasis basis_q, basis_x;
98 CeedVector q_data, x_coord, jac_data = NULL;
99 CeedOperator op_ifunction_vol = NULL, op_rhs_vol = NULL, op_ijacobian_vol = NULL;
100
101 PetscFunctionBeginUser;
102 PetscCall(DMGetDimension(dm, &dim));
103 PetscCall(DMGetCoordinateNumComps(dm, &num_comp_x));
104
105 CeedElemRestriction elem_restr_diff_flux = NULL;
106 CeedVector div_diff_flux_ceed = NULL;
107 CeedBasis basis_diff_flux = NULL;
108 CeedEvalMode eval_mode_diff_flux = -1;
109 { // Create bases and element restrictions
110 PetscInt height = 0, dm_field = 0;
111 DM dm_coord;
112
113 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
114 PetscCall(DMPlexCeedBasisCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, dm_field, &basis_q));
115
116 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, 0, &elem_restr_q));
117 if (num_comps_jac_data) {
118 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, num_comps_jac_data, &elem_restr_jd_i));
119 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_jd_i, &jac_data, NULL));
120 }
121
122 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q, &honee->q_ceed, NULL));
123 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q, &honee->q_dot_ceed, NULL));
124 PetscCallCeed(ceed, CeedElemRestrictionCreateVector(elem_restr_q, &honee->g_ceed, NULL));
125 PetscCall(DMPlexCeedCoordinateCreateField(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, &elem_restr_x, &basis_x, &x_coord));
126
127 PetscCall(QDataGet(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size_vol));
128 }
129
130 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) {
131 PetscCheck(honee->diff_flux_proj, honee->comm, PETSC_ERR_ARG_WRONGSTATE,
132 "Divergence of diffusive flux projection requested but object not created");
133 PetscCall(DivDiffFluxProjectionGetOperatorFieldData(honee->diff_flux_proj, &elem_restr_diff_flux, &basis_diff_flux, &div_diff_flux_ceed,
134 &eval_mode_diff_flux));
135 }
136
137 { // -- Create QFunction for ICs
138 CeedBasis basis_xc;
139 CeedQFunction qf_ics;
140 CeedOperator op_ics;
141
142 PetscCallCeed(ceed, CeedBasisCreateProjection(basis_x, basis_q, &basis_xc));
143 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->ics.qf_func_ptr, problem->ics.qf_loc, &qf_ics));
144 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ics, problem->ics.qfctx));
145 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ics, 0));
146 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ics, "x", num_comp_x, CEED_EVAL_INTERP));
147 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ics, "dx", num_comp_x * dim, CEED_EVAL_GRAD));
148 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ics, "q0", num_comp_q, CEED_EVAL_NONE));
149
150 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics));
151 PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "x", elem_restr_x, basis_xc, CEED_VECTOR_ACTIVE));
152 PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "dx", elem_restr_x, basis_xc, CEED_VECTOR_ACTIVE));
153 PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "q0", elem_restr_q, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
154 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_ics, "evaluation time", &honee->phys->ics_time_label));
155 PetscCall(OperatorApplyContextCreate(NULL, dm, honee->ceed, op_ics, x_coord, NULL, NULL, honee->Q_loc, &honee->op_ics_ctx));
156
157 PetscCallCeed(ceed, CeedBasisDestroy(&basis_xc));
158 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ics));
159 PetscCallCeed(ceed, CeedOperatorDestroy(&op_ics));
160 }
161
162 if (problem->apply_vol_rhs.qf_func_ptr) {
163 CeedQFunction qf_rhs_vol;
164
165 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs.qf_func_ptr, problem->apply_vol_rhs.qf_loc, &qf_rhs_vol));
166 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_vol, problem->apply_vol_rhs.qfctx));
167 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_rhs_vol, 0));
168 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP));
169 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD));
170 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
171 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP));
172 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
173 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "div F_diff", honee->diff_flux_proj->num_diff_flux_comps, eval_mode_diff_flux));
174 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_vol, "v", num_comp_q, CEED_EVAL_INTERP));
175 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
176
177 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_vol, NULL, NULL, &op_rhs_vol));
178 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
179 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
180 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
181 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "x", elem_restr_x, basis_x, x_coord));
182 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
183 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "div F_diff", elem_restr_diff_flux, basis_diff_flux, div_diff_flux_ceed));
184 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
185 PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
186
187 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_vol));
188 }
189
190 if (problem->apply_vol_ifunction.qf_func_ptr) {
191 CeedQFunction qf_ifunction_vol;
192
193 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction.qf_func_ptr, problem->apply_vol_ifunction.qf_loc,
194 &qf_ifunction_vol));
195 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ifunction_vol, problem->apply_vol_ifunction.qfctx));
196 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ifunction_vol, 0));
197 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "q", num_comp_q, CEED_EVAL_INTERP));
198 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD));
199 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "q dot", num_comp_q, CEED_EVAL_INTERP));
200 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
201 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "x", num_comp_x, CEED_EVAL_INTERP));
202 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
203 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "div F_diff", honee->diff_flux_proj->num_diff_flux_comps, eval_mode_diff_flux));
204 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "v", num_comp_q, CEED_EVAL_INTERP));
205 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
206 if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "jac_data", num_comps_jac_data, CEED_EVAL_NONE));
207
208 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ifunction_vol, NULL, NULL, &op_ifunction_vol));
209 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
210 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
211 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "q dot", elem_restr_q, basis_q, honee->q_dot_ceed));
212 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
213 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "x", elem_restr_x, basis_x, x_coord));
214 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
215 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "div F_diff", elem_restr_diff_flux, basis_diff_flux, div_diff_flux_ceed));
216 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
217 PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
218 if (num_comps_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "jac_data", elem_restr_jd_i, CEED_BASIS_NONE, jac_data));
219
220 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ifunction_vol));
221 }
222
223 if (problem->apply_vol_ijacobian.qf_func_ptr) {
224 CeedQFunction qf_ijacobian_vol;
225
226 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ijacobian.qf_func_ptr, problem->apply_vol_ijacobian.qf_loc,
227 &qf_ijacobian_vol));
228 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ijacobian_vol, problem->apply_vol_ijacobian.qfctx));
229 PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ijacobian_vol, 0));
230 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "dq", num_comp_q, CEED_EVAL_INTERP));
231 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "Grad_dq", num_comp_q * dim, CEED_EVAL_GRAD));
232 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
233 if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "jac_data", num_comps_jac_data, CEED_EVAL_NONE));
234 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ijacobian_vol, "v", num_comp_q, CEED_EVAL_INTERP));
235 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ijacobian_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
236
237 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ijacobian_vol, NULL, NULL, &op_ijacobian_vol));
238 PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
239 PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "Grad_dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
240 PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
241 if (num_comps_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "jac_data", elem_restr_jd_i, CEED_BASIS_NONE, jac_data));
242 PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
243 PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
244
245 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ijacobian_vol));
246 }
247
248 // -- Create and apply CEED Composite Operator for the entire domain
249 if (!honee->phys->implicit) { // RHS
250 CeedOperator op_rhs;
251
252 PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_rhs));
253 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_rhs, op_rhs_vol));
254 for (PetscCount b = 0; b < problem->num_bc_defs; b++) PetscCall(BCDefinitionAddOperators(problem->bc_defs[b], op_rhs, NULL));
255
256 PetscCall(OperatorApplyContextCreate(dm, dm, ceed, op_rhs, honee->q_ceed, honee->g_ceed, honee->Q_loc, NULL, &honee->op_rhs_ctx));
257
258 // ----- Get Context Labels for Operator
259 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_rhs, "solution time", &honee->phys->solution_time_label));
260 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_rhs, "timestep size", &honee->phys->timestep_size_label));
261
262 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
263 PetscCall(CreateKSPMass(honee, problem));
264 PetscCheck(app_ctx->sgs_model_type == SGS_MODEL_NONE, honee->comm, PETSC_ERR_SUP, "SGS modeling not implemented for explicit timestepping");
265 } else { // IFunction
266 CeedOperator op_ijacobian = NULL;
267
268 // Create Composite Operaters
269 PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &honee->op_ifunction));
270 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(honee->op_ifunction, op_ifunction_vol));
271 if (op_ijacobian_vol) {
272 PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_ijacobian));
273 PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ijacobian, op_ijacobian_vol));
274 }
275 for (PetscCount b = 0; b < problem->num_bc_defs; b++) PetscCall(BCDefinitionAddOperators(problem->bc_defs[b], honee->op_ifunction, op_ijacobian));
276
277 // ----- Get Context Labels for Operator
278 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(honee->op_ifunction, "solution time", &honee->phys->solution_time_label));
279 PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(honee->op_ifunction, "timestep size", &honee->phys->timestep_size_label));
280
281 if (op_ijacobian) {
282 PetscCall(MatCreateCeed(honee->dm, honee->dm, op_ijacobian, NULL, &honee->mat_ijacobian));
283 PetscCall(MatCeedSetLocalVectors(honee->mat_ijacobian, honee->Q_dot_loc, NULL));
284 PetscCallCeed(ceed, CeedOperatorDestroy(&op_ijacobian));
285 }
286 if (app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) PetscCall(SgsDDSetup(ceed, honee, problem));
287 }
288
289 if (problem->use_strong_bc_ceed) PetscCall(SetupStrongBC_Ceed(ceed, dm, honee, problem));
290 if (app_ctx->sgs_train_enable) PetscCall(SGS_DD_TrainingSetup(ceed, honee));
291 if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) PetscCall(DivDiffFluxProjectionSetup(honee, honee->diff_flux_proj));
292
293 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
294 PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
295 PetscCallCeed(ceed, CeedVectorDestroy(&x_coord));
296 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
297 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_x));
298 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
299 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jd_i));
300 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
301 PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
302 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
303 PetscCallCeed(ceed, CeedBasisDestroy(&basis_x));
304 PetscCallCeed(ceed, CeedVectorDestroy(&div_diff_flux_ceed));
305 PetscCallCeed(ceed, CeedOperatorDestroy(&op_ijacobian_vol));
306 PetscCallCeed(ceed, CeedOperatorDestroy(&op_ifunction_vol));
307 PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs_vol));
308 PetscFunctionReturn(PETSC_SUCCESS);
309 }
310