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