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 /// @file
9 /// libCEED setup for solid mechanics example using PETSc
10
11 #include "../include/setup-libceed.h"
12
13 #include <ceed.h>
14 #include <petscdmplex.h>
15 #include <petscmat.h>
16
17 #include "../include/structs.h"
18 #include "../include/utils.h"
19 #include "../qfunctions/constant-force.h" // Constant forcing function
20 #include "../qfunctions/manufactured-force.h" // Manufactured solution forcing
21 #include "../qfunctions/traction-boundary.h" // Traction boundaries
22
23 // -----------------------------------------------------------------------------
24 // Problem options
25 // -----------------------------------------------------------------------------
26 // Forcing function data
27 forcingData forcing_options[3] = {
28 [FORCE_NONE] = {.setup_forcing = NULL, .setup_forcing_loc = NULL },
29 [FORCE_CONST] = {.setup_forcing = SetupConstantForce, .setup_forcing_loc = SetupConstantForce_loc},
30 [FORCE_MMS] = {.setup_forcing = SetupMMSForce, .setup_forcing_loc = SetupMMSForce_loc }
31 };
32
33 // -----------------------------------------------------------------------------
34 // libCEED Functions
35 // -----------------------------------------------------------------------------
36 // Destroy libCEED objects
CeedDataDestroy(CeedInt level,CeedData data)37 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data) {
38 PetscFunctionBegin;
39
40 // Vectors
41 CeedVectorDestroy(&data->x_ceed);
42 CeedVectorDestroy(&data->y_ceed);
43 CeedVectorDestroy(&data->geo_data);
44 for (CeedInt i = 0; i < SOLIDS_MAX_NUMBER_FIELDS; i++) CeedVectorDestroy(&data->stored_fields[i]);
45 CeedVectorDestroy(&data->geo_data_diagnostic);
46 CeedVectorDestroy(&data->true_soln);
47 // Restrictions
48 CeedElemRestrictionDestroy(&data->elem_restr_x);
49 CeedElemRestrictionDestroy(&data->elem_restr_u);
50 CeedElemRestrictionDestroy(&data->elem_restr_geo_data_i);
51 for (CeedInt i = 0; i < SOLIDS_MAX_NUMBER_FIELDS; i++) CeedElemRestrictionDestroy(&data->elem_restr_stored_fields_i[i]);
52 CeedElemRestrictionDestroy(&data->elem_restr_energy);
53 CeedElemRestrictionDestroy(&data->elem_restr_diagnostic);
54 CeedElemRestrictionDestroy(&data->elem_restr_geo_data_diagnostic_i);
55 // Bases
56 CeedBasisDestroy(&data->basis_x);
57 CeedBasisDestroy(&data->basis_u);
58 CeedBasisDestroy(&data->basis_energy);
59 CeedBasisDestroy(&data->basis_diagnostic);
60 // QFunctions
61 CeedQFunctionDestroy(&data->qf_residual);
62 CeedQFunctionDestroy(&data->qf_jacobian);
63 CeedQFunctionDestroy(&data->qf_energy);
64 CeedQFunctionDestroy(&data->qf_diagnostic);
65 // Operators
66 CeedOperatorDestroy(&data->op_residual);
67 CeedOperatorDestroy(&data->op_jacobian);
68 CeedOperatorDestroy(&data->op_energy);
69 CeedOperatorDestroy(&data->op_diagnostic);
70 // Restriction and Prolongation data
71 CeedBasisDestroy(&data->basis_c_to_f);
72 CeedOperatorDestroy(&data->op_prolong);
73 CeedOperatorDestroy(&data->op_restrict);
74
75 PetscCall(PetscFree(data));
76
77 PetscFunctionReturn(PETSC_SUCCESS);
78 };
79
80 // Utility function to create local CEED restriction from DMPlex
CreateRestrictionFromPlex(Ceed ceed,DM dm,CeedInt height,DMLabel domain_label,CeedInt value,CeedElemRestriction * elem_restr)81 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) {
82 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets_petsc;
83 CeedInt *elem_restr_offsets_ceed;
84
85 PetscFunctionBeginUser;
86 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_petsc));
87
88 PetscCall(IntArrayPetscToCeed(num_elem * elem_size, &elem_restr_offsets_petsc, &elem_restr_offsets_ceed));
89 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets_ceed, elem_restr);
90 PetscCall(PetscFree(elem_restr_offsets_ceed));
91
92 PetscFunctionReturn(PETSC_SUCCESS);
93 };
94
95 // Utility function to get Ceed Restriction for each domain
GetRestrictionForDomain(Ceed ceed,DM dm,CeedInt height,DMLabel domain_label,PetscInt value,CeedInt Q,CeedInt q_data_size,CeedElemRestriction * elem_restr_q,CeedElemRestriction * elem_restr_x,CeedElemRestriction * elem_restr_qd_i)96 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, PetscInt value, CeedInt Q, CeedInt q_data_size,
97 CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, CeedElemRestriction *elem_restr_qd_i) {
98 DM dm_coord;
99 CeedInt num_local_elem, Q_dim;
100 PetscInt dim;
101
102 PetscFunctionBeginUser;
103
104 PetscCall(DMGetDimension(dm, &dim));
105 dim -= height;
106 Q_dim = CeedIntPow(Q, dim);
107 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
108 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
109 if (elem_restr_q) {
110 PetscCall(CreateRestrictionFromPlex(ceed, dm, height, domain_label, value, elem_restr_q));
111 }
112 if (elem_restr_x) {
113 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, height, domain_label, value, elem_restr_x));
114 }
115 if (elem_restr_qd_i) {
116 CeedElemRestrictionGetNumElements(*elem_restr_q, &num_local_elem);
117 CeedElemRestrictionCreateStrided(ceed, num_local_elem, Q_dim, q_data_size, q_data_size * num_local_elem * Q_dim, CEED_STRIDES_BACKEND,
118 elem_restr_qd_i);
119 }
120
121 PetscFunctionReturn(PETSC_SUCCESS);
122 };
123
124 // Set up libCEED on the fine grid for a given degree
SetupLibceedFineLevel(DM dm,DM dm_energy,DM dm_diagnostic,Ceed ceed,AppCtx app_ctx,CeedQFunctionContext phys_ctx,ProblemData problem_data,PetscInt fine_level,PetscInt num_comp_u,PetscInt U_g_size,PetscInt U_loc_size,CeedVector force_ceed,CeedVector neumann_ceed,CeedData * data)125 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dm_energy, DM dm_diagnostic, Ceed ceed, AppCtx app_ctx, CeedQFunctionContext phys_ctx,
126 ProblemData problem_data, PetscInt fine_level, PetscInt num_comp_u, PetscInt U_g_size, PetscInt U_loc_size,
127 CeedVector force_ceed, CeedVector neumann_ceed, CeedData *data) {
128 CeedInt P = app_ctx->level_degrees[fine_level] + 1;
129 CeedInt Q = app_ctx->level_degrees[fine_level] + 1 + app_ctx->q_extra;
130 CeedInt num_comp_x, num_comp_e = 1, num_comp_d = 5;
131 CeedInt num_qpts;
132 CeedInt q_data_size = problem_data.q_data_size;
133 forcingType forcing_choice = app_ctx->forcing_choice;
134 DM dm_coord;
135 Vec coords;
136 PetscInt c_start, c_end, num_elem, dim;
137 const PetscScalar *coordArray;
138 CeedVector x_coord;
139 CeedQFunction qf_setup_geo, qf_residual, qf_jacobian, qf_energy, qf_diagnostic;
140 CeedOperator op_setup_geo, op_residual, op_jacobian, op_energy, op_diagnostic;
141
142 PetscFunctionBeginUser;
143
144 // ---------------------------------------------------------------------------
145 // libCEED bases
146 // ---------------------------------------------------------------------------
147 PetscCall(DMGetDimension(dm, &dim));
148 num_comp_x = dim;
149 // -- Coordinate basis
150 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, problem_data.quadrature_mode, &data[fine_level]->basis_x);
151 // -- Solution basis
152 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, problem_data.quadrature_mode, &data[fine_level]->basis_u);
153 // -- Energy basis
154 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_e, P, Q, problem_data.quadrature_mode, &data[fine_level]->basis_energy);
155 // -- Diagnostic output basis
156 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, P, CEED_GAUSS_LOBATTO, &data[fine_level]->basis_diagnostic);
157
158 // ---------------------------------------------------------------------------
159 // libCEED restrictions
160 // ---------------------------------------------------------------------------
161 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
162 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
163
164 // -- Coordinate restriction
165 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &(data[fine_level]->elem_restr_x)));
166 // -- Solution restriction
167 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &data[fine_level]->elem_restr_u));
168 // -- Energy restriction
169 PetscCall(CreateRestrictionFromPlex(ceed, dm_energy, 0, 0, 0, &data[fine_level]->elem_restr_energy));
170 // -- Diagnostic data restriction
171 PetscCall(CreateRestrictionFromPlex(ceed, dm_diagnostic, 0, 0, 0, &data[fine_level]->elem_restr_diagnostic));
172
173 // -- Stored data at quadrature points
174 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end));
175 num_elem = c_end - c_start;
176 CeedBasisGetNumQuadraturePoints(data[fine_level]->basis_u, &num_qpts);
177 // ---- Geometric data restriction, residual and Jacobian operators
178 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, num_elem * num_qpts * q_data_size, CEED_STRIDES_BACKEND,
179 &data[fine_level]->elem_restr_geo_data_i);
180 // ---- Stored field restrictions
181 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
182 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, problem_data.field_sizes[i], num_elem * num_qpts * problem_data.field_sizes[i],
183 CEED_STRIDES_BACKEND, &data[fine_level]->elem_restr_stored_fields_i[i]);
184 }
185 // ---- Geometric data restriction, diagnostic operator
186 CeedElemRestrictionCreateStrided(ceed, num_elem, P * P * P, q_data_size, num_elem * P * P * P * q_data_size, CEED_STRIDES_BACKEND,
187 &data[fine_level]->elem_restr_geo_data_diagnostic_i);
188
189 // ---------------------------------------------------------------------------
190 // Element coordinates
191 // ---------------------------------------------------------------------------
192 PetscCall(DMGetCoordinatesLocal(dm, &coords));
193 PetscCall(VecGetArrayRead(coords, &coordArray));
194
195 CeedElemRestrictionCreateVector(data[fine_level]->elem_restr_x, &x_coord, NULL);
196 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coordArray);
197 PetscCall(VecRestoreArrayRead(coords, &coordArray));
198
199 // ---------------------------------------------------------------------------
200 // Persistent libCEED vectors
201 // ---------------------------------------------------------------------------
202 // -- Operator action variables
203 CeedVectorCreate(ceed, U_loc_size, &data[fine_level]->x_ceed);
204 CeedVectorCreate(ceed, U_loc_size, &data[fine_level]->y_ceed);
205 // -- Geometric data vector
206 CeedVectorCreate(ceed, num_elem * num_qpts * q_data_size, &data[fine_level]->geo_data);
207 // -- Stored field vectors
208 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
209 CeedVectorCreate(ceed, num_elem * num_qpts * problem_data.field_sizes[i], &data[fine_level]->stored_fields[i]);
210 }
211 // -- Collocated geometric data vector
212 CeedVectorCreate(ceed, num_elem * P * P * P * q_data_size, &data[fine_level]->geo_data_diagnostic);
213
214 // ---------------------------------------------------------------------------
215 // Geometric factor computation
216 // ---------------------------------------------------------------------------
217 // Create the QFunction and Operator that computes the quadrature data geo_data returns dXdx_i,j and w * det.
218 // ---------------------------------------------------------------------------
219 // -- QFunction
220 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, problem_data.setup_geo_loc, &qf_setup_geo);
221 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * dim, CEED_EVAL_GRAD);
222 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT);
223 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE);
224 // -- Operator
225 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_geo);
226 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, data[fine_level]->basis_x, CEED_VECTOR_ACTIVE);
227 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, data[fine_level]->basis_x, CEED_VECTOR_NONE);
228 CeedOperatorSetField(op_setup_geo, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
229 // -- Compute the quadrature data
230 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data, CEED_REQUEST_IMMEDIATE);
231 // -- Cleanup
232 CeedQFunctionDestroy(&qf_setup_geo);
233 CeedOperatorDestroy(&op_setup_geo);
234
235 // ---------------------------------------------------------------------------
236 // Local residual evaluator
237 // ---------------------------------------------------------------------------
238 // Create the QFunction and Operator that computes the residual of the non-linear PDE.
239 // ---------------------------------------------------------------------------
240 // -- QFunction
241 CeedQFunctionCreateInterior(ceed, 1, problem_data.residual, problem_data.residual_loc, &qf_residual);
242 CeedQFunctionAddInput(qf_residual, "du", num_comp_u * dim, CEED_EVAL_GRAD);
243 CeedQFunctionAddInput(qf_residual, "qdata", q_data_size, CEED_EVAL_NONE);
244 CeedQFunctionAddOutput(qf_residual, "dv", num_comp_u * dim, CEED_EVAL_GRAD);
245 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
246 CeedQFunctionAddOutput(qf_residual, problem_data.field_names[i], problem_data.field_sizes[i], CEED_EVAL_NONE);
247 }
248 CeedQFunctionSetContext(qf_residual, phys_ctx);
249 // -- Operator
250 CeedOperatorCreate(ceed, qf_residual, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_residual);
251 CeedOperatorSetField(op_residual, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
252 CeedOperatorSetField(op_residual, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_NONE, data[fine_level]->geo_data);
253 CeedOperatorSetField(op_residual, "dv", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
254 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
255 CeedOperatorSetField(op_residual, problem_data.field_names[i], data[fine_level]->elem_restr_stored_fields_i[i], CEED_BASIS_NONE,
256 data[fine_level]->stored_fields[i]);
257 }
258 // -- Save libCEED data
259 data[fine_level]->qf_residual = qf_residual;
260 data[fine_level]->op_residual = op_residual;
261
262 // ---------------------------------------------------------------------------
263 // Jacobian evaluator
264 // ---------------------------------------------------------------------------
265 // Create the QFunction and Operator that computes the action of the Jacobian for each linear solve.
266 // ---------------------------------------------------------------------------
267 // -- QFunction
268 CeedQFunctionCreateInterior(ceed, 1, problem_data.jacobian, problem_data.jacobian_loc, &qf_jacobian);
269 CeedQFunctionAddInput(qf_jacobian, "delta du", num_comp_u * dim, CEED_EVAL_GRAD);
270 CeedQFunctionAddInput(qf_jacobian, "qdata", q_data_size, CEED_EVAL_NONE);
271 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
272 CeedQFunctionAddInput(qf_jacobian, problem_data.field_names[i], problem_data.field_sizes[i], CEED_EVAL_NONE);
273 }
274 CeedQFunctionAddOutput(qf_jacobian, "delta dv", num_comp_u * dim, CEED_EVAL_GRAD);
275 CeedQFunctionSetContext(qf_jacobian, phys_ctx);
276 // -- Operator
277 CeedOperatorCreate(ceed, qf_jacobian, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_jacobian);
278 CeedOperatorSetField(op_jacobian, "delta du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
279 CeedOperatorSetField(op_jacobian, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_NONE, data[fine_level]->geo_data);
280 CeedOperatorSetField(op_jacobian, "delta dv", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
281 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) {
282 CeedOperatorSetField(op_jacobian, problem_data.field_names[i], data[fine_level]->elem_restr_stored_fields_i[i], CEED_BASIS_NONE,
283 data[fine_level]->stored_fields[i]);
284 }
285 // -- Save libCEED data
286 data[fine_level]->qf_jacobian = qf_jacobian;
287 data[fine_level]->op_jacobian = op_jacobian;
288
289 // ---------------------------------------------------------------------------
290 // Traction boundary conditions, if needed
291 // ---------------------------------------------------------------------------
292 if (app_ctx->bc_traction_count > 0) {
293 // -- Setup
294 DMLabel domain_label;
295 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label));
296 PetscCall(DMGetDimension(dm, &dim));
297
298 // -- Basis
299 CeedInt height = 1;
300 CeedBasis basis_x_face, basis_u_face;
301 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_x, 2, Q, problem_data.quadrature_mode, &basis_x_face);
302 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_u, P, Q, problem_data.quadrature_mode, &basis_u_face);
303 // -- QFunction
304 CeedQFunction qf_traction;
305 CeedQFunctionContext traction_ctx;
306 CeedQFunctionCreateInterior(ceed, 1, SetupTractionBCs, SetupTractionBCs_loc, &qf_traction);
307 CeedQFunctionContextCreate(ceed, &traction_ctx);
308 CeedQFunctionSetContext(qf_traction, traction_ctx);
309 CeedQFunctionAddInput(qf_traction, "dx", num_comp_x * (num_comp_x - height), CEED_EVAL_GRAD);
310 CeedQFunctionAddInput(qf_traction, "weight", 1, CEED_EVAL_WEIGHT);
311 CeedQFunctionAddOutput(qf_traction, "v", num_comp_u, CEED_EVAL_INTERP);
312
313 // -- Compute contribution on each boundary face
314 for (CeedInt i = 0; i < app_ctx->bc_traction_count; i++) {
315 CeedElemRestriction elem_restr_x_face, elem_restr_u_face;
316 CeedOperator op_traction;
317 CeedQFunctionContextSetData(traction_ctx, CEED_MEM_HOST, CEED_USE_POINTER, 3 * sizeof(CeedScalar), app_ctx->bc_traction_vector[i]);
318 // Setup restriction
319 PetscCall(GetRestrictionForDomain(ceed, dm, 1, domain_label, app_ctx->bc_traction_faces[i], Q, 0, &elem_restr_u_face, &elem_restr_x_face,
320 NULL));
321 // ---- Create boundary Operator
322 CeedOperatorCreate(ceed, qf_traction, NULL, NULL, &op_traction);
323 CeedOperatorSetField(op_traction, "dx", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE);
324 CeedOperatorSetField(op_traction, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE);
325 CeedOperatorSetField(op_traction, "v", elem_restr_u_face, basis_u_face, CEED_VECTOR_ACTIVE);
326 // ---- Compute traction on face
327 CeedOperatorApplyAdd(op_traction, x_coord, neumann_ceed, CEED_REQUEST_IMMEDIATE);
328 // ---- Cleanup
329 CeedElemRestrictionDestroy(&elem_restr_x_face);
330 CeedElemRestrictionDestroy(&elem_restr_u_face);
331 CeedOperatorDestroy(&op_traction);
332 }
333 // -- Cleanup
334 CeedBasisDestroy(&basis_x_face);
335 CeedBasisDestroy(&basis_u_face);
336 CeedQFunctionDestroy(&qf_traction);
337 CeedQFunctionContextDestroy(&traction_ctx);
338 }
339
340 // ---------------------------------------------------------------------------
341 // Forcing term, if needed
342 // ---------------------------------------------------------------------------
343 // Create the QFunction and Operator that computes the forcing term (RHS) for the non-linear PDE.
344 // ---------------------------------------------------------------------------
345 if (forcing_choice != FORCE_NONE) {
346 CeedQFunction qf_setup_force;
347 CeedOperator op_setup_force;
348
349 // -- QFunction
350 CeedQFunctionCreateInterior(ceed, 1, forcing_options[forcing_choice].setup_forcing, forcing_options[forcing_choice].setup_forcing_loc,
351 &qf_setup_force);
352 CeedQFunctionAddInput(qf_setup_force, "x", num_comp_x, CEED_EVAL_INTERP);
353 CeedQFunctionAddInput(qf_setup_force, "qdata", q_data_size, CEED_EVAL_NONE);
354 CeedQFunctionAddOutput(qf_setup_force, "force", num_comp_u, CEED_EVAL_INTERP);
355 if (forcing_choice == FORCE_MMS) {
356 CeedQFunctionSetContext(qf_setup_force, phys_ctx);
357 } else {
358 CeedQFunctionContext ctxForcing;
359 CeedQFunctionContextCreate(ceed, &ctxForcing);
360 CeedQFunctionContextSetData(ctxForcing, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*app_ctx->forcing_vector), app_ctx->forcing_vector);
361 CeedQFunctionSetContext(qf_setup_force, ctxForcing);
362 CeedQFunctionContextDestroy(&ctxForcing);
363 }
364 // -- Operator
365 CeedOperatorCreate(ceed, qf_setup_force, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_force);
366 CeedOperatorSetField(op_setup_force, "x", data[fine_level]->elem_restr_x, data[fine_level]->basis_x, CEED_VECTOR_ACTIVE);
367 CeedOperatorSetField(op_setup_force, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_NONE, data[fine_level]->geo_data);
368 CeedOperatorSetField(op_setup_force, "force", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
369 // -- Compute forcing term
370 CeedOperatorApply(op_setup_force, x_coord, force_ceed, CEED_REQUEST_IMMEDIATE);
371 // -- Cleanup
372 CeedQFunctionDestroy(&qf_setup_force);
373 CeedOperatorDestroy(&op_setup_force);
374 }
375
376 // ---------------------------------------------------------------------------
377 // True solution, for MMS
378 // ---------------------------------------------------------------------------
379 // Create the QFunction and Operator that computes the true solution at the mesh nodes for validation with the manufactured solution.
380 // ---------------------------------------------------------------------------
381 if (problem_data.true_soln) {
382 CeedScalar *true_array;
383 const CeedScalar *mult_array;
384 CeedVector mult_vec;
385 CeedBasis basis_x_true;
386 CeedQFunction qf_true;
387 CeedOperator op_true;
388
389 // -- Solution vector
390 CeedVectorCreate(ceed, U_loc_size, &(data[fine_level]->true_soln));
391 // -- Basis
392 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, CEED_GAUSS_LOBATTO, &basis_x_true);
393 // QFunction
394 CeedQFunctionCreateInterior(ceed, 1, problem_data.true_soln, problem_data.true_soln_loc, &qf_true);
395 CeedQFunctionAddInput(qf_true, "x", num_comp_x, CEED_EVAL_INTERP);
396 CeedQFunctionAddOutput(qf_true, "true solution", num_comp_u, CEED_EVAL_NONE);
397 // Operator
398 CeedOperatorCreate(ceed, qf_true, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_true);
399 CeedOperatorSetField(op_true, "x", data[fine_level]->elem_restr_x, basis_x_true, CEED_VECTOR_ACTIVE);
400 CeedOperatorSetField(op_true, "true solution", data[fine_level]->elem_restr_u, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
401 // -- Compute true solution
402 CeedOperatorApply(op_true, x_coord, data[fine_level]->true_soln, CEED_REQUEST_IMMEDIATE);
403 // -- Multiplicity calculation
404 CeedElemRestrictionCreateVector(data[fine_level]->elem_restr_u, &mult_vec, NULL);
405 CeedVectorSetValue(mult_vec, 0.);
406 CeedElemRestrictionGetMultiplicity(data[fine_level]->elem_restr_u, mult_vec);
407 // -- Multiplicity correction
408 CeedVectorGetArray(data[fine_level]->true_soln, CEED_MEM_HOST, &true_array);
409 CeedVectorGetArrayRead(mult_vec, CEED_MEM_HOST, &mult_array);
410 for (CeedInt i = 0; i < U_loc_size; i++) true_array[i] /= mult_array[i];
411 CeedVectorRestoreArray(data[fine_level]->true_soln, &true_array);
412 CeedVectorRestoreArrayRead(mult_vec, &mult_array);
413 // -- Cleanup
414 CeedVectorDestroy(&mult_vec);
415 CeedBasisDestroy(&basis_x_true);
416 CeedQFunctionDestroy(&qf_true);
417 CeedOperatorDestroy(&op_true);
418 }
419
420 // ---------------------------------------------------------------------------
421 // Local energy computation
422 // ---------------------------------------------------------------------------
423 // Create the QFunction and Operator that computes the strain energy
424 // ---------------------------------------------------------------------------
425 // -- QFunction
426 CeedQFunctionCreateInterior(ceed, 1, problem_data.energy, problem_data.energy_loc, &qf_energy);
427 CeedQFunctionAddInput(qf_energy, "du", num_comp_u * dim, CEED_EVAL_GRAD);
428 CeedQFunctionAddInput(qf_energy, "qdata", q_data_size, CEED_EVAL_NONE);
429 CeedQFunctionAddOutput(qf_energy, "energy", num_comp_e, CEED_EVAL_INTERP);
430 CeedQFunctionSetContext(qf_energy, phys_ctx);
431 // -- Operator
432 CeedOperatorCreate(ceed, qf_energy, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_energy);
433 CeedOperatorSetField(op_energy, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE);
434 CeedOperatorSetField(op_energy, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_NONE, data[fine_level]->geo_data);
435 CeedOperatorSetField(op_energy, "energy", data[fine_level]->elem_restr_energy, data[fine_level]->basis_energy, CEED_VECTOR_ACTIVE);
436 // -- Save libCEED data
437 data[fine_level]->qf_energy = qf_energy;
438 data[fine_level]->op_energy = op_energy;
439
440 // ---------------------------------------------------------------------------
441 // Diagnostic value computation
442 // ---------------------------------------------------------------------------
443 // Create the QFunction and Operator that computes nodal diagnostic quantities
444 // ---------------------------------------------------------------------------
445 // Geometric factors
446 // -- Coordinate basis
447 CeedBasis basis_x;
448 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS_LOBATTO, &basis_x);
449 // -- QFunction
450 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, problem_data.setup_geo_loc, &qf_setup_geo);
451 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * dim, CEED_EVAL_GRAD);
452 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT);
453 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE);
454 // -- Operator
455 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_geo);
456 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, basis_x, CEED_VECTOR_ACTIVE);
457 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
458 CeedOperatorSetField(op_setup_geo, "qdata", data[fine_level]->elem_restr_geo_data_diagnostic_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
459 // -- Compute the quadrature data
460 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data_diagnostic, CEED_REQUEST_IMMEDIATE);
461 // -- Cleanup
462 CeedBasisDestroy(&basis_x);
463 CeedQFunctionDestroy(&qf_setup_geo);
464 CeedOperatorDestroy(&op_setup_geo);
465
466 // Diagnostic quantities
467 // -- QFunction
468 CeedQFunctionCreateInterior(ceed, 1, problem_data.diagnostic, problem_data.diagnostic_loc, &qf_diagnostic);
469 CeedQFunctionAddInput(qf_diagnostic, "u", num_comp_u, CEED_EVAL_INTERP);
470 CeedQFunctionAddInput(qf_diagnostic, "du", num_comp_u * dim, CEED_EVAL_GRAD);
471 CeedQFunctionAddInput(qf_diagnostic, "qdata", q_data_size, CEED_EVAL_NONE);
472 CeedQFunctionAddOutput(qf_diagnostic, "diagnostic values", num_comp_u + num_comp_d, CEED_EVAL_NONE);
473 CeedQFunctionSetContext(qf_diagnostic, phys_ctx);
474 // -- Operator
475 CeedOperatorCreate(ceed, qf_diagnostic, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diagnostic);
476 CeedOperatorSetField(op_diagnostic, "u", data[fine_level]->elem_restr_u, data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE);
477 CeedOperatorSetField(op_diagnostic, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE);
478 CeedOperatorSetField(op_diagnostic, "qdata", data[fine_level]->elem_restr_geo_data_diagnostic_i, CEED_BASIS_NONE,
479 data[fine_level]->geo_data_diagnostic);
480 CeedOperatorSetField(op_diagnostic, "diagnostic values", data[fine_level]->elem_restr_diagnostic, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
481 // -- Save libCEED data
482 data[fine_level]->qf_diagnostic = qf_diagnostic;
483 data[fine_level]->op_diagnostic = op_diagnostic;
484
485 // ---------------------------------------------------------------------------
486 // Cleanup
487 // ---------------------------------------------------------------------------
488 CeedVectorDestroy(&x_coord);
489
490 PetscFunctionReturn(PETSC_SUCCESS);
491 };
492
493 // Set up libCEED multigrid level for a given degree
494 // Prolongation and Restriction are between level and level+1
SetupLibceedLevel(DM dm,Ceed ceed,AppCtx app_ctx,ProblemData problem_data,PetscInt level,PetscInt num_comp_u,PetscInt U_g_size,PetscInt U_loc_size,CeedVector fine_mult,CeedData * data)495 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData problem_data, PetscInt level, PetscInt num_comp_u, PetscInt U_g_size,
496 PetscInt U_loc_size, CeedVector fine_mult, CeedData *data) {
497 CeedInt fine_level = app_ctx->num_levels - 1;
498 CeedInt P = app_ctx->level_degrees[level] + 1;
499 CeedInt Q = app_ctx->level_degrees[fine_level] + 1 + app_ctx->q_extra;
500 PetscInt dim;
501 CeedOperator op_jacobian, op_prolong, op_restrict;
502
503 PetscFunctionBeginUser;
504
505 PetscCall(DMGetDimension(dm, &dim));
506
507 // ---------------------------------------------------------------------------
508 // libCEED restrictions
509 // ---------------------------------------------------------------------------
510 // -- Solution restriction
511 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &data[level]->elem_restr_u));
512
513 // ---------------------------------------------------------------------------
514 // libCEED bases
515 // ---------------------------------------------------------------------------
516 // -- Solution basis
517 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, problem_data.quadrature_mode, &data[level]->basis_u);
518
519 // ---------------------------------------------------------------------------
520 // Persistent libCEED vectors
521 // ---------------------------------------------------------------------------
522 CeedVectorCreate(ceed, U_loc_size, &data[level]->x_ceed);
523 CeedVectorCreate(ceed, U_loc_size, &data[level]->y_ceed);
524
525 // ---------------------------------------------------------------------------
526 // Coarse Grid, Prolongation, and Restriction Operators
527 // ---------------------------------------------------------------------------
528 // Create the Operators that compute the prolongation and restriction between the p-multigrid levels and the coarse grid eval.
529 // ---------------------------------------------------------------------------
530 CeedOperatorMultigridLevelCreate(data[level + 1]->op_jacobian, fine_mult, data[level]->elem_restr_u, data[level]->basis_u, &op_jacobian,
531 &op_prolong, &op_restrict);
532
533 // -- Save libCEED data
534 data[level]->op_jacobian = op_jacobian;
535 data[level + 1]->op_prolong = op_prolong;
536 data[level + 1]->op_restrict = op_restrict;
537
538 PetscFunctionReturn(PETSC_SUCCESS);
539 };
540