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