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 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 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 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 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 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