xref: /honee/src/setuplibceed.c (revision 77efe802949112fcee8af9d4702ab5d884d2f004)
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
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
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 
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             num_comps_jac_data = problem->num_comps_jac_data, q_data_size_vol;
96   CeedElemRestriction elem_restr_jd_i    = NULL, elem_restr_qd, elem_restr_q;
97   CeedBasis           basis_q;
98   CeedVector          jac_data         = NULL, q_data;
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 
126     PetscCall(DMPlexCeedCoordinateGetField(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, height, &honee->elem_restr_x, &honee->basis_x,
127                                            &honee->x_coord));
128 
129     PetscCall(QDataGet(ceed, dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, &q_data,
130                        &q_data_size_vol));
131   }
132 
133   if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) {
134     PetscCheck(honee->diff_flux_proj, honee->comm, PETSC_ERR_ARG_WRONGSTATE,
135                "Divergence of diffusive flux projection requested but object not created");
136     PetscCall(DivDiffFluxProjectionGetOperatorFieldData(honee->diff_flux_proj, &elem_restr_diff_flux, &basis_diff_flux, &div_diff_flux_ceed,
137                                                         &eval_mode_diff_flux));
138   }
139 
140   {  // -- Create QFunction for ICs
141     CeedBasis     basis_xc;
142     CeedQFunction qf_ics;
143     CeedOperator  op_ics;
144 
145     PetscCallCeed(ceed, CeedBasisCreateProjection(honee->basis_x, basis_q, &basis_xc));
146     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->ics.qf_func_ptr, problem->ics.qf_loc, &qf_ics));
147     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ics, problem->ics.qfctx));
148     PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ics, 0));
149     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ics, "x", num_comp_x, CEED_EVAL_INTERP));
150     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ics, "dx", num_comp_x * dim, CEED_EVAL_GRAD));
151     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ics, "q0", num_comp_q, CEED_EVAL_NONE));
152 
153     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ics, NULL, NULL, &op_ics));
154     PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "x", honee->elem_restr_x, basis_xc, CEED_VECTOR_ACTIVE));
155     PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "dx", honee->elem_restr_x, basis_xc, CEED_VECTOR_ACTIVE));
156     PetscCallCeed(ceed, CeedOperatorSetField(op_ics, "q0", elem_restr_q, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
157     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_ics, "evaluation time", &honee->phys->ics_time_label));
158     PetscCall(OperatorApplyContextCreate(NULL, dm, honee->ceed, op_ics, honee->x_coord, NULL, NULL, honee->Q_loc, &honee->op_ics_ctx));
159 
160     PetscCallCeed(ceed, CeedBasisDestroy(&basis_xc));
161     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ics));
162     PetscCallCeed(ceed, CeedOperatorDestroy(&op_ics));
163   }
164 
165   if (problem->apply_vol_rhs.qf_func_ptr) {
166     CeedQFunction qf_rhs_vol;
167 
168     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs.qf_func_ptr, problem->apply_vol_rhs.qf_loc, &qf_rhs_vol));
169     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_rhs_vol, problem->apply_vol_rhs.qfctx));
170     PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_rhs_vol, 0));
171     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "q", num_comp_q, CEED_EVAL_INTERP));
172     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD));
173     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
174     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "x", num_comp_x, CEED_EVAL_INTERP));
175     if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
176       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_rhs_vol, "div F_diff", honee->diff_flux_proj->num_diff_flux_comps, eval_mode_diff_flux));
177     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_vol, "v", num_comp_q, CEED_EVAL_INTERP));
178     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_rhs_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
179 
180     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_rhs_vol, NULL, NULL, &op_rhs_vol));
181     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
182     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
183     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
184     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "x", honee->elem_restr_x, honee->basis_x, honee->x_coord));
185     if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
186       PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "div F_diff", elem_restr_diff_flux, basis_diff_flux, div_diff_flux_ceed));
187     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
188     PetscCallCeed(ceed, CeedOperatorSetField(op_rhs_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
189 
190     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_rhs_vol));
191   }
192 
193   if (problem->apply_vol_ifunction.qf_func_ptr) {
194     CeedQFunction qf_ifunction_vol;
195 
196     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ifunction.qf_func_ptr, problem->apply_vol_ifunction.qf_loc,
197                                                     &qf_ifunction_vol));
198     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ifunction_vol, problem->apply_vol_ifunction.qfctx));
199     PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ifunction_vol, 0));
200     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "q", num_comp_q, CEED_EVAL_INTERP));
201     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD));
202     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "q dot", num_comp_q, CEED_EVAL_INTERP));
203     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
204     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "x", num_comp_x, CEED_EVAL_INTERP));
205     if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
206       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ifunction_vol, "div F_diff", honee->diff_flux_proj->num_diff_flux_comps, eval_mode_diff_flux));
207     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "v", num_comp_q, CEED_EVAL_INTERP));
208     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
209     if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ifunction_vol, "jac_data", num_comps_jac_data, CEED_EVAL_NONE));
210 
211     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ifunction_vol, NULL, NULL, &op_ifunction_vol));
212     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
213     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
214     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "q dot", elem_restr_q, basis_q, honee->q_dot_ceed));
215     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
216     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "x", honee->elem_restr_x, honee->basis_x, honee->x_coord));
217     if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE)
218       PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "div F_diff", elem_restr_diff_flux, basis_diff_flux, div_diff_flux_ceed));
219     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
220     PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
221     if (num_comps_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(op_ifunction_vol, "jac_data", elem_restr_jd_i, CEED_BASIS_NONE, jac_data));
222 
223     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ifunction_vol));
224   }
225 
226   if (problem->apply_vol_ijacobian.qf_func_ptr) {
227     CeedQFunction qf_ijacobian_vol;
228 
229     PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_ijacobian.qf_func_ptr, problem->apply_vol_ijacobian.qf_loc,
230                                                     &qf_ijacobian_vol));
231     PetscCallCeed(ceed, CeedQFunctionSetContext(qf_ijacobian_vol, problem->apply_vol_ijacobian.qfctx));
232     PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_ijacobian_vol, 0));
233     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "dq", num_comp_q, CEED_EVAL_INTERP));
234     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "Grad_dq", num_comp_q * dim, CEED_EVAL_GRAD));
235     PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "qdata", q_data_size_vol, CEED_EVAL_NONE));
236     if (num_comps_jac_data) PetscCallCeed(ceed, CeedQFunctionAddInput(qf_ijacobian_vol, "jac_data", num_comps_jac_data, CEED_EVAL_NONE));
237     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ijacobian_vol, "v", num_comp_q, CEED_EVAL_INTERP));
238     PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_ijacobian_vol, "Grad_v", num_comp_q * dim, CEED_EVAL_GRAD));
239 
240     PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_ijacobian_vol, NULL, NULL, &op_ijacobian_vol));
241     PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
242     PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "Grad_dq", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
243     PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "qdata", elem_restr_qd, CEED_BASIS_NONE, q_data));
244     if (num_comps_jac_data) PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "jac_data", elem_restr_jd_i, CEED_BASIS_NONE, jac_data));
245     PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
246     PetscCallCeed(ceed, CeedOperatorSetField(op_ijacobian_vol, "Grad_v", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
247 
248     PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_ijacobian_vol));
249   }
250 
251   // -- Create and apply CEED Composite Operator for the entire domain
252   if (!honee->phys->implicit) {  // RHS
253     CeedOperator op_rhs;
254 
255     PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_rhs));
256     PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_rhs, op_rhs_vol));
257     for (PetscCount b = 0; b < problem->num_bc_defs; b++) PetscCall(BCDefinitionAddOperators(problem->bc_defs[b], op_rhs, NULL));
258 
259     PetscCall(OperatorApplyContextCreate(dm, dm, ceed, op_rhs, honee->q_ceed, honee->g_ceed, honee->Q_loc, NULL, &honee->op_rhs_ctx));
260 
261     // ----- Get Context Labels for Operator
262     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_rhs, "solution time", &honee->phys->solution_time_label));
263     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(op_rhs, "timestep size", &honee->phys->timestep_size_label));
264 
265     PetscCallCeed(ceed, CeedOperatorDestroy(&op_rhs));
266     PetscCall(CreateKSPMass(honee, problem));
267     PetscCheck(app_ctx->sgs_model_type == SGS_MODEL_NONE, honee->comm, PETSC_ERR_SUP, "SGS modeling not implemented for explicit timestepping");
268   } else {  // IFunction
269     CeedOperator op_ijacobian = NULL;
270 
271     // Create Composite Operaters
272     PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &honee->op_ifunction));
273     PetscCallCeed(ceed, CeedOperatorCompositeAddSub(honee->op_ifunction, op_ifunction_vol));
274     if (op_ijacobian_vol) {
275       PetscCallCeed(ceed, CeedOperatorCreateComposite(ceed, &op_ijacobian));
276       PetscCallCeed(ceed, CeedOperatorCompositeAddSub(op_ijacobian, op_ijacobian_vol));
277     }
278     for (PetscCount b = 0; b < problem->num_bc_defs; b++) PetscCall(BCDefinitionAddOperators(problem->bc_defs[b], honee->op_ifunction, op_ijacobian));
279 
280     // ----- Get Context Labels for Operator
281     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(honee->op_ifunction, "solution time", &honee->phys->solution_time_label));
282     PetscCallCeed(ceed, CeedOperatorGetContextFieldLabel(honee->op_ifunction, "timestep size", &honee->phys->timestep_size_label));
283 
284     if (op_ijacobian) {
285       PetscCall(MatCreateCeed(honee->dm, honee->dm, op_ijacobian, NULL, &honee->mat_ijacobian));
286       PetscCall(MatCeedSetLocalVectors(honee->mat_ijacobian, honee->Q_dot_loc, NULL));
287       PetscCallCeed(ceed, CeedOperatorDestroy(&op_ijacobian));
288     }
289     if (app_ctx->sgs_model_type == SGS_MODEL_DATA_DRIVEN) PetscCall(SgsDDSetup(ceed, honee, problem));
290   }
291 
292   if (problem->use_strong_bc_ceed) PetscCall(SetupStrongBC_Ceed(ceed, dm, honee, problem));
293   if (app_ctx->sgs_train_enable) PetscCall(SGS_DD_TrainingSetup(ceed, honee));
294   if (app_ctx->divFdiffproj_method != DIV_DIFF_FLUX_PROJ_NONE) PetscCall(DivDiffFluxProjectionSetup(honee, honee->diff_flux_proj));
295 
296   PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
297   PetscCallCeed(ceed, CeedVectorDestroy(&jac_data));
298   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
299   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
300   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_jd_i));
301   PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_diff_flux));
302   PetscCallCeed(ceed, CeedBasisDestroy(&basis_diff_flux));
303   PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
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