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