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 - essential BC dofs are encoded in closure indices as -(i+1) 87 PetscInt Involute(PetscInt i) { 88 return i >= 0 ? i : -(i + 1); 89 }; 90 91 // Utility function to create local CEED restriction from DMPlex 92 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 93 CeedInt height, DMLabel domain_label, CeedInt value, 94 CeedElemRestriction *elem_restr) { 95 96 PetscSection section; 97 PetscInt p, num_elem, num_dof, *restr_indices, elem_offset, num_fields, dim, 98 depth; 99 DMLabel depth_label; 100 IS depth_is, iter_is; 101 Vec U_loc; 102 const PetscInt *iter_indices; 103 PetscErrorCode ierr; 104 105 PetscFunctionBeginUser; 106 107 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 108 dim -= height; 109 ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr); 110 ierr = PetscSectionGetNumFields(section, &num_fields); CHKERRQ(ierr); 111 PetscInt num_comp[num_fields], field_offsets[num_fields+1]; 112 field_offsets[0] = 0; 113 for (PetscInt f = 0; f < num_fields; f++) { 114 ierr = PetscSectionGetFieldComponents(section, f, &num_comp[f]); CHKERRQ(ierr); 115 field_offsets[f+1] = field_offsets[f] + num_comp[f]; 116 } 117 118 ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); 119 ierr = DMPlexGetDepthLabel(dm, &depth_label); CHKERRQ(ierr); 120 ierr = DMLabelGetStratumIS(depth_label, depth - height, &depth_is); 121 CHKERRQ(ierr); 122 if (domain_label) { 123 IS domain_is; 124 ierr = DMLabelGetStratumIS(domain_label, value, &domain_is); CHKERRQ(ierr); 125 if (domain_is) { // domainIS is non-empty 126 ierr = ISIntersect(depth_is, domain_is, &iter_is); CHKERRQ(ierr); 127 ierr = ISDestroy(&domain_is); CHKERRQ(ierr); 128 } else { // domainIS is NULL (empty) 129 iter_is = NULL; 130 } 131 ierr = ISDestroy(&depth_is); CHKERRQ(ierr); 132 } else { 133 iter_is = depth_is; 134 } 135 if (iter_is) { 136 ierr = ISGetLocalSize(iter_is, &num_elem); CHKERRQ(ierr); 137 ierr = ISGetIndices(iter_is, &iter_indices); CHKERRQ(ierr); 138 } else { 139 num_elem = 0; 140 iter_indices = NULL; 141 } 142 ierr = PetscMalloc1(num_elem*PetscPowInt(P, dim), &restr_indices); 143 CHKERRQ(ierr); 144 for (p = 0, elem_offset = 0; p < num_elem; p++) { 145 PetscInt c = iter_indices[p]; 146 PetscInt num_indices, *indices, num_nodes; 147 ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, 148 &num_indices, &indices, NULL, NULL); 149 CHKERRQ(ierr); 150 bool flip = false; 151 if (height > 0) { 152 PetscInt num_cells, num_faces, start = -1; 153 const PetscInt *orients, *faces, *cells; 154 ierr = DMPlexGetSupport(dm, c, &cells); CHKERRQ(ierr); 155 ierr = DMPlexGetSupportSize(dm, c, &num_cells); CHKERRQ(ierr); 156 if (num_cells != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 157 "Expected one cell in support of exterior face, but got %D cells", 158 num_cells); 159 ierr = DMPlexGetCone(dm, cells[0], &faces); CHKERRQ(ierr); 160 ierr = DMPlexGetConeSize(dm, cells[0], &num_faces); CHKERRQ(ierr); 161 for (PetscInt i=0; i<num_faces; i++) {if (faces[i] == c) start = i;} 162 if (start < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, 163 "Could not find face %D in cone of its support", 164 c); 165 ierr = DMPlexGetConeOrientation(dm, cells[0], &orients); CHKERRQ(ierr); 166 if (orients[start] < 0) flip = true; 167 } 168 if (num_indices % field_offsets[num_fields]) SETERRQ1(PETSC_COMM_SELF, 169 PETSC_ERR_ARG_INCOMP, "Number of closure indices not compatible with Cell %D", 170 c); 171 num_nodes = num_indices / field_offsets[num_fields]; 172 for (PetscInt i = 0; i < num_nodes; i++) { 173 PetscInt ii = i; 174 if (flip) { 175 if (P == num_nodes) ii = num_nodes - 1 - i; 176 else if (P*P == num_nodes) { 177 PetscInt row = i / P, col = i % P; 178 ii = row + col * P; 179 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, 180 "No support for flipping point with %D nodes != P (%D) or P^2", 181 num_nodes, P); 182 } 183 // Check that indices are blocked by node and thus can be coalesced as a single field with 184 // field_offsets[num_fields] = sum(num_comp) components. 185 for (PetscInt f = 0; f < num_fields; f++) { 186 for (PetscInt j = 0; j < num_comp[f]; j++) { 187 if (Involute(indices[field_offsets[f]*num_nodes + ii*num_comp[f] + j]) 188 != Involute(indices[ii*num_comp[0]]) + field_offsets[f] + j) 189 SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 190 "Cell %D closure indices not interlaced for node %D field %D component %D", 191 c, ii, f, j); 192 } 193 } 194 // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. 195 PetscInt loc = Involute(indices[ii*num_comp[0]]); 196 restr_indices[elem_offset++] = loc; 197 } 198 ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, 199 &num_indices, &indices, NULL, NULL); 200 CHKERRQ(ierr); 201 } 202 if (elem_offset != num_elem*PetscPowInt(P, dim)) 203 SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, 204 "ElemRestriction of size (%D,%D) initialized %D nodes", num_elem, 205 PetscPowInt(P, dim),elem_offset); 206 if (iter_is) { 207 ierr = ISRestoreIndices(iter_is, &iter_indices); CHKERRQ(ierr); 208 } 209 ierr = ISDestroy(&iter_is); CHKERRQ(ierr); 210 211 ierr = DMGetLocalVector(dm, &U_loc); CHKERRQ(ierr); 212 ierr = VecGetLocalSize(U_loc, &num_dof); CHKERRQ(ierr); 213 ierr = DMRestoreLocalVector(dm, &U_loc); CHKERRQ(ierr); 214 CeedElemRestrictionCreate(ceed, num_elem, PetscPowInt(P, dim), 215 field_offsets[num_fields], 216 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices, 217 elem_restr); 218 ierr = PetscFree(restr_indices); CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 }; 221 222 // Utility function to get Ceed Restriction for each domain 223 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 224 DMLabel domain_label, PetscInt value, 225 CeedInt P, CeedInt Q, CeedInt q_data_size, 226 CeedElemRestriction *elem_restr_q, 227 CeedElemRestriction *elem_restr_x, 228 CeedElemRestriction *elem_restr_qd_i) { 229 230 DM dm_coord; 231 CeedInt dim, num_local_elem; 232 CeedInt Q_dim; 233 PetscErrorCode ierr; 234 235 PetscFunctionBeginUser; 236 237 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 238 dim -= height; 239 Q_dim = CeedIntPow(Q, dim); 240 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 241 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 242 CHKERRQ(ierr); 243 if (elem_restr_q) { 244 ierr = CreateRestrictionFromPlex(ceed, dm, P, height, domain_label, value, 245 elem_restr_q); CHKERRQ(ierr); 246 } 247 if (elem_restr_x) { 248 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 2, height, domain_label, 249 value, elem_restr_x); CHKERRQ(ierr); 250 } 251 if (elem_restr_qd_i) { 252 CeedElemRestrictionGetNumElements(*elem_restr_q, &num_local_elem); 253 CeedElemRestrictionCreateStrided(ceed, num_local_elem, Q_dim, 254 q_data_size, q_data_size*num_local_elem*Q_dim, 255 CEED_STRIDES_BACKEND, elem_restr_qd_i); 256 } 257 258 PetscFunctionReturn(0); 259 }; 260 261 // Set up libCEED on the fine grid for a given degree 262 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dm_energy, DM dm_diagnostic, 263 Ceed ceed, AppCtx app_ctx, 264 CeedQFunctionContext phys_ctx, 265 ProblemData problem_data, 266 PetscInt fine_level, PetscInt num_comp_u, 267 PetscInt U_g_size, PetscInt U_loc_size, 268 CeedVector force_ceed, 269 CeedVector neumann_ceed, CeedData *data) { 270 int ierr; 271 CeedInt P = app_ctx->level_degrees[fine_level] + 1; 272 CeedInt Q = app_ctx->level_degrees[fine_level] + 1 + app_ctx->q_extra; 273 CeedInt dim, num_comp_x, num_comp_e = 1, num_comp_d = 5; 274 CeedInt num_qpts; 275 CeedInt geo_data_size = problem_data.geo_data_size; 276 forcingType forcing_choice = app_ctx->forcing_choice; 277 DM dm_coord; 278 Vec coords; 279 PetscInt c_start, c_end, num_elem; 280 const PetscScalar *coordArray; 281 CeedVector x_coord; 282 CeedQFunction qf_setup_geo, qf_residual, qf_jacobian, qf_energy, qf_diagnostic; 283 CeedOperator op_setup_geo, op_residual, op_jacobian, op_energy, op_diagnostic; 284 285 PetscFunctionBeginUser; 286 287 // --------------------------------------------------------------------------- 288 // libCEED bases 289 // --------------------------------------------------------------------------- 290 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 291 num_comp_x = dim; 292 // -- Coordinate basis 293 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, 294 problem_data.quadrature_mode, 295 &data[fine_level]->basis_x); 296 // -- Solution basis 297 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, 298 problem_data.quadrature_mode, 299 &data[fine_level]->basis_u); 300 // -- Energy basis 301 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_e, P, Q, 302 problem_data.quadrature_mode, 303 &data[fine_level]->basis_energy); 304 // -- Diagnostic output basis 305 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, P, CEED_GAUSS_LOBATTO, 306 &data[fine_level]->basis_diagnostic); 307 308 // --------------------------------------------------------------------------- 309 // libCEED restrictions 310 // --------------------------------------------------------------------------- 311 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 312 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 313 CHKERRQ(ierr); 314 315 // -- Coordinate restriction 316 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 2, 0, 0, 0, 317 &(data[fine_level]->elem_restr_x)); 318 CHKERRQ(ierr); 319 // -- Solution restriction 320 ierr = CreateRestrictionFromPlex(ceed, dm, P, 0, 0, 0, 321 &data[fine_level]->elem_restr_u); 322 CHKERRQ(ierr); 323 // -- Energy restriction 324 ierr = CreateRestrictionFromPlex(ceed, dm_energy, P, 0, 0, 0, 325 &data[fine_level]->elem_restr_energy); 326 CHKERRQ(ierr); 327 // -- Diagnostic data restriction 328 ierr = CreateRestrictionFromPlex(ceed, dm_diagnostic, P, 0, 0, 0, 329 &data[fine_level]->elem_restr_diagnostic); 330 CHKERRQ(ierr); 331 332 // -- Stored data at quadrature points 333 ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); CHKERRQ(ierr); 334 num_elem = c_end - c_start; 335 CeedBasisGetNumQuadraturePoints(data[fine_level]->basis_u, &num_qpts); 336 // ---- Geometric data restriction, residual and Jacobian operators 337 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, geo_data_size, 338 num_elem*num_qpts*geo_data_size, 339 CEED_STRIDES_BACKEND, 340 &data[fine_level]->elem_restr_geo_data_i); 341 // ---- Stored field restrictions 342 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 343 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, 344 problem_data.field_sizes[i], 345 num_elem*num_qpts*problem_data.field_sizes[i], 346 CEED_STRIDES_BACKEND, 347 &data[fine_level]->elem_restr_stored_fields_i[i]); 348 } 349 // ---- Geometric data restriction, diagnostic operator 350 CeedElemRestrictionCreateStrided(ceed, num_elem, P*P*P, geo_data_size, 351 num_elem*P*P*P*geo_data_size, 352 CEED_STRIDES_BACKEND, 353 &data[fine_level]->elem_restr_geo_data_diagnostic_i); 354 355 // --------------------------------------------------------------------------- 356 // Element coordinates 357 // --------------------------------------------------------------------------- 358 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 359 ierr = VecGetArrayRead(coords, &coordArray); CHKERRQ(ierr); 360 361 CeedElemRestrictionCreateVector(data[fine_level]->elem_restr_x, &x_coord, NULL); 362 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 363 (PetscScalar *)coordArray); 364 ierr = VecRestoreArrayRead(coords, &coordArray); CHKERRQ(ierr); 365 366 // --------------------------------------------------------------------------- 367 // Persistent libCEED vectors 368 // --------------------------------------------------------------------------- 369 // -- Operator action variables 370 CeedVectorCreate(ceed, U_loc_size, &data[fine_level]->x_ceed); 371 CeedVectorCreate(ceed, U_loc_size, &data[fine_level]->y_ceed); 372 // -- Geometric data vector 373 CeedVectorCreate(ceed, num_elem*num_qpts*geo_data_size, 374 &data[fine_level]->geo_data); 375 // -- Stored field vectors 376 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 377 CeedVectorCreate(ceed, num_elem*num_qpts*problem_data.field_sizes[i], 378 &data[fine_level]->stored_fields[i]); 379 } 380 // -- Collocated geometric data vector 381 CeedVectorCreate(ceed, num_elem*P*P*P*geo_data_size, 382 &data[fine_level]->geo_data_diagnostic); 383 384 // --------------------------------------------------------------------------- 385 // Geometric factor computation 386 // --------------------------------------------------------------------------- 387 // Create the QFunction and Operator that computes the quadrature data 388 // geo_data returns dXdx_i,j and w * det. 389 // --------------------------------------------------------------------------- 390 // -- QFunction 391 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, 392 problem_data.setup_geo_loc, &qf_setup_geo); 393 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*dim, CEED_EVAL_GRAD); 394 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 395 CeedQFunctionAddOutput(qf_setup_geo, "geo_data", geo_data_size, CEED_EVAL_NONE); 396 // -- Operator 397 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, 398 CEED_QFUNCTION_NONE, &op_setup_geo); 399 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, 400 data[fine_level]->basis_x, CEED_VECTOR_ACTIVE); 401 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, 402 data[fine_level]->basis_x, CEED_VECTOR_NONE); 403 CeedOperatorSetField(op_setup_geo, "geo_data", 404 data[fine_level]->elem_restr_geo_data_i, 405 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 406 // -- Compute the quadrature data 407 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data, 408 CEED_REQUEST_IMMEDIATE); 409 // -- Cleanup 410 CeedQFunctionDestroy(&qf_setup_geo); 411 CeedOperatorDestroy(&op_setup_geo); 412 413 // --------------------------------------------------------------------------- 414 // Local residual evaluator 415 // --------------------------------------------------------------------------- 416 // Create the QFunction and Operator that computes the residual of the 417 // non-linear PDE. 418 // --------------------------------------------------------------------------- 419 // -- QFunction 420 CeedQFunctionCreateInterior(ceed, 1, problem_data.residual, 421 problem_data.residual_loc, &qf_residual); 422 CeedQFunctionAddInput(qf_residual, "du", num_comp_u*dim, CEED_EVAL_GRAD); 423 CeedQFunctionAddInput(qf_residual, "geo_data", geo_data_size, CEED_EVAL_NONE); 424 CeedQFunctionAddOutput(qf_residual, "dv", num_comp_u*dim, CEED_EVAL_GRAD); 425 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 426 CeedQFunctionAddOutput(qf_residual, problem_data.field_names[i], 427 problem_data.field_sizes[i], CEED_EVAL_NONE); 428 } 429 CeedQFunctionSetContext(qf_residual, phys_ctx); 430 // -- Operator 431 CeedOperatorCreate(ceed, qf_residual, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 432 &op_residual); 433 CeedOperatorSetField(op_residual, "du", data[fine_level]->elem_restr_u, 434 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 435 CeedOperatorSetField(op_residual, "geo_data", 436 data[fine_level]->elem_restr_geo_data_i, 437 CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 438 CeedOperatorSetField(op_residual, "dv", data[fine_level]->elem_restr_u, 439 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 440 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 441 CeedOperatorSetField(op_residual, problem_data.field_names[i], 442 data[fine_level]->elem_restr_stored_fields_i[i], 443 CEED_BASIS_COLLOCATED, 444 data[fine_level]->stored_fields[i]); 445 } 446 // -- Save libCEED data 447 data[fine_level]->qf_residual = qf_residual; 448 data[fine_level]->op_residual = op_residual; 449 450 // --------------------------------------------------------------------------- 451 // Jacobian evaluator 452 // --------------------------------------------------------------------------- 453 // Create the QFunction and Operator that computes the action of the 454 // Jacobian for each linear solve. 455 // --------------------------------------------------------------------------- 456 // -- QFunction 457 CeedQFunctionCreateInterior(ceed, 1, problem_data.jacobian, 458 problem_data.jacobian_loc, &qf_jacobian); 459 CeedQFunctionAddInput(qf_jacobian, "deltadu", num_comp_u*dim, CEED_EVAL_GRAD); 460 CeedQFunctionAddInput(qf_jacobian, "geo_data", geo_data_size, CEED_EVAL_NONE); 461 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 462 CeedQFunctionAddInput(qf_jacobian, problem_data.field_names[i], 463 problem_data.field_sizes[i], CEED_EVAL_NONE); 464 } 465 CeedQFunctionAddOutput(qf_jacobian, "deltadv", num_comp_u*dim, CEED_EVAL_GRAD); 466 CeedQFunctionSetContext(qf_jacobian, phys_ctx); 467 // -- Operator 468 CeedOperatorCreate(ceed, qf_jacobian, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 469 &op_jacobian); 470 CeedOperatorSetField(op_jacobian, "deltadu", data[fine_level]->elem_restr_u, 471 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 472 CeedOperatorSetField(op_jacobian, "geo_data", 473 data[fine_level]->elem_restr_geo_data_i, 474 CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 475 CeedOperatorSetField(op_jacobian, "deltadv", data[fine_level]->elem_restr_u, 476 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 477 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 478 CeedOperatorSetField(op_jacobian, problem_data.field_names[i], 479 data[fine_level]->elem_restr_stored_fields_i[i], 480 CEED_BASIS_COLLOCATED, 481 data[fine_level]->stored_fields[i]); 482 } 483 // -- Save libCEED data 484 data[fine_level]->qf_jacobian = qf_jacobian; 485 data[fine_level]->op_jacobian = op_jacobian; 486 487 // --------------------------------------------------------------------------- 488 // Traction boundary conditions, if needed 489 // --------------------------------------------------------------------------- 490 if (app_ctx->bc_traction_count > 0) { 491 // -- Setup 492 DMLabel domain_label; 493 ierr = DMGetLabel(dm, "Face Sets", &domain_label); CHKERRQ(ierr); 494 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 495 496 // -- Basis 497 CeedInt height = 1; 498 CeedBasis basis_x_face, basis_u_face; 499 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_x, 2, Q, 500 problem_data.quadrature_mode, &basis_x_face); 501 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_u, P, Q, 502 problem_data.quadrature_mode, &basis_u_face); 503 // -- QFunction 504 CeedQFunction qf_traction; 505 CeedQFunctionContext traction_ctx; 506 CeedQFunctionCreateInterior(ceed, 1, SetupTractionBCs, SetupTractionBCs_loc, 507 &qf_traction); 508 CeedQFunctionContextCreate(ceed, &traction_ctx); 509 CeedQFunctionSetContext(qf_traction, traction_ctx); 510 CeedQFunctionAddInput(qf_traction, "dx", num_comp_x*(num_comp_x - height), 511 CEED_EVAL_GRAD); 512 CeedQFunctionAddInput(qf_traction, "weight", 1, CEED_EVAL_WEIGHT); 513 CeedQFunctionAddOutput(qf_traction, "v", num_comp_u, CEED_EVAL_INTERP); 514 515 // -- Compute contribution on each boundary face 516 for (CeedInt i = 0; i < app_ctx->bc_traction_count; i++) { 517 CeedElemRestriction elem_restr_x_face, elem_restr_u_face; 518 CeedOperator op_traction; 519 CeedQFunctionContextSetData(traction_ctx, CEED_MEM_HOST, CEED_USE_POINTER, 520 3 * sizeof(CeedScalar), 521 app_ctx->bc_traction_vector[i]); 522 // Setup restriction 523 ierr = GetRestrictionForDomain(ceed, dm, height, domain_label, 524 app_ctx->bc_traction_faces[i], P, Q, 525 0, &elem_restr_u_face, &elem_restr_x_face, NULL); 526 CHKERRQ(ierr); 527 // ---- Create boundary Operator 528 CeedOperatorCreate(ceed, qf_traction, NULL, NULL, &op_traction); 529 CeedOperatorSetField(op_traction, "dx", elem_restr_x_face, basis_x_face, 530 CEED_VECTOR_ACTIVE); 531 CeedOperatorSetField(op_traction, "weight", CEED_ELEMRESTRICTION_NONE, 532 basis_x_face, CEED_VECTOR_NONE); 533 CeedOperatorSetField(op_traction, "v", elem_restr_u_face, 534 basis_u_face, CEED_VECTOR_ACTIVE); 535 // ---- Compute traction on face 536 CeedOperatorApplyAdd(op_traction, x_coord, neumann_ceed, 537 CEED_REQUEST_IMMEDIATE); 538 // ---- Cleanup 539 CeedElemRestrictionDestroy(&elem_restr_x_face); 540 CeedElemRestrictionDestroy(&elem_restr_u_face); 541 CeedOperatorDestroy(&op_traction); 542 } 543 // -- Cleanup 544 CeedBasisDestroy(&basis_x_face); 545 CeedBasisDestroy(&basis_u_face); 546 CeedQFunctionDestroy(&qf_traction); 547 CeedQFunctionContextDestroy(&traction_ctx); 548 } 549 550 // --------------------------------------------------------------------------- 551 // Forcing term, if needed 552 // --------------------------------------------------------------------------- 553 // Create the QFunction and Operator that computes the forcing term (RHS) 554 // for the non-linear PDE. 555 // --------------------------------------------------------------------------- 556 if (forcing_choice != FORCE_NONE) { 557 CeedQFunction qf_setup_force; 558 CeedOperator op_setup_force; 559 560 // -- QFunction 561 CeedQFunctionCreateInterior(ceed, 1, 562 forcing_options[forcing_choice].setup_forcing, 563 forcing_options[forcing_choice].setup_forcing_loc, 564 &qf_setup_force); 565 CeedQFunctionAddInput(qf_setup_force, "x", num_comp_x, CEED_EVAL_INTERP); 566 CeedQFunctionAddInput(qf_setup_force, "geo_data", geo_data_size, 567 CEED_EVAL_NONE); 568 CeedQFunctionAddOutput(qf_setup_force, "force", num_comp_u, CEED_EVAL_INTERP); 569 if (forcing_choice == FORCE_MMS) { 570 CeedQFunctionSetContext(qf_setup_force, phys_ctx); 571 } else { 572 CeedQFunctionContext ctxForcing; 573 CeedQFunctionContextCreate(ceed, &ctxForcing); 574 CeedQFunctionContextSetData(ctxForcing, CEED_MEM_HOST, CEED_USE_POINTER, 575 sizeof(*app_ctx->forcing_vector), 576 app_ctx->forcing_vector); 577 CeedQFunctionSetContext(qf_setup_force, ctxForcing); 578 CeedQFunctionContextDestroy(&ctxForcing); 579 } 580 // -- Operator 581 CeedOperatorCreate(ceed, qf_setup_force, CEED_QFUNCTION_NONE, 582 CEED_QFUNCTION_NONE, &op_setup_force); 583 CeedOperatorSetField(op_setup_force, "x", data[fine_level]->elem_restr_x, 584 data[fine_level]->basis_x, CEED_VECTOR_ACTIVE); 585 CeedOperatorSetField(op_setup_force, "geo_data", 586 data[fine_level]->elem_restr_geo_data_i, 587 CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 588 CeedOperatorSetField(op_setup_force, "force", data[fine_level]->elem_restr_u, 589 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 590 // -- Compute forcing term 591 CeedOperatorApply(op_setup_force, x_coord, force_ceed, CEED_REQUEST_IMMEDIATE); 592 // -- Cleanup 593 CeedQFunctionDestroy(&qf_setup_force); 594 CeedOperatorDestroy(&op_setup_force); 595 } 596 597 // --------------------------------------------------------------------------- 598 // True solution, for MMS 599 // --------------------------------------------------------------------------- 600 // Create the QFunction and Operator that computes the true solution at 601 // the mesh nodes for validation with the manufactured solution. 602 // --------------------------------------------------------------------------- 603 if (problem_data.true_soln) { 604 CeedScalar *true_array; 605 const CeedScalar *mult_array; 606 CeedVector mult_vec; 607 CeedBasis basis_x_true; 608 CeedQFunction qf_true; 609 CeedOperator op_true; 610 611 // -- Solution vector 612 CeedVectorCreate(ceed, U_loc_size, &(data[fine_level]->true_soln)); 613 // -- Basis 614 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, CEED_GAUSS_LOBATTO, 615 &basis_x_true); 616 // QFunction 617 CeedQFunctionCreateInterior(ceed, 1, problem_data.true_soln, 618 problem_data.true_soln_loc, &qf_true); 619 CeedQFunctionAddInput(qf_true, "x", num_comp_x, CEED_EVAL_INTERP); 620 CeedQFunctionAddOutput(qf_true, "true_soln", num_comp_u, CEED_EVAL_NONE); 621 // Operator 622 CeedOperatorCreate(ceed, qf_true, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 623 &op_true); 624 CeedOperatorSetField(op_true, "x", data[fine_level]->elem_restr_x, basis_x_true, 625 CEED_VECTOR_ACTIVE); 626 CeedOperatorSetField(op_true, "true_soln", data[fine_level]->elem_restr_u, 627 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 628 // -- Compute true solution 629 CeedOperatorApply(op_true, x_coord, data[fine_level]->true_soln, 630 CEED_REQUEST_IMMEDIATE); 631 // -- Multiplicity calculation 632 CeedElemRestrictionCreateVector(data[fine_level]->elem_restr_u, &mult_vec, 633 NULL); 634 CeedVectorSetValue(mult_vec, 0.); 635 CeedElemRestrictionGetMultiplicity(data[fine_level]->elem_restr_u, mult_vec); 636 // -- Multiplicity correction 637 CeedVectorGetArray(data[fine_level]->true_soln, CEED_MEM_HOST, &true_array); 638 CeedVectorGetArrayRead(mult_vec, CEED_MEM_HOST, &mult_array); 639 for (CeedInt i = 0; i < U_loc_size; i++) 640 true_array[i] /= mult_array[i]; 641 CeedVectorRestoreArray(data[fine_level]->true_soln, &true_array); 642 CeedVectorRestoreArrayRead(mult_vec, &mult_array); 643 // -- Cleanup 644 CeedVectorDestroy(&mult_vec); 645 CeedBasisDestroy(&basis_x_true); 646 CeedQFunctionDestroy(&qf_true); 647 CeedOperatorDestroy(&op_true); 648 } 649 650 // --------------------------------------------------------------------------- 651 // Local energy computation 652 // --------------------------------------------------------------------------- 653 // Create the QFunction and Operator that computes the strain energy 654 // --------------------------------------------------------------------------- 655 // -- QFunction 656 CeedQFunctionCreateInterior(ceed, 1, problem_data.energy, 657 problem_data.energy_loc, &qf_energy); 658 CeedQFunctionAddInput(qf_energy, "du", num_comp_u*dim, CEED_EVAL_GRAD); 659 CeedQFunctionAddInput(qf_energy, "geo_data", geo_data_size, CEED_EVAL_NONE); 660 CeedQFunctionAddOutput(qf_energy, "energy", num_comp_e, CEED_EVAL_INTERP); 661 CeedQFunctionSetContext(qf_energy, phys_ctx); 662 // -- Operator 663 CeedOperatorCreate(ceed, qf_energy, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 664 &op_energy); 665 CeedOperatorSetField(op_energy, "du", data[fine_level]->elem_restr_u, 666 data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 667 CeedOperatorSetField(op_energy, "geo_data", 668 data[fine_level]->elem_restr_geo_data_i, 669 CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 670 CeedOperatorSetField(op_energy, "energy", data[fine_level]->elem_restr_energy, 671 data[fine_level]->basis_energy, CEED_VECTOR_ACTIVE); 672 // -- Save libCEED data 673 data[fine_level]->qf_energy = qf_energy; 674 data[fine_level]->op_energy = op_energy; 675 676 // --------------------------------------------------------------------------- 677 // Diagnostic value computation 678 // --------------------------------------------------------------------------- 679 // Create the QFunction and Operator that computes nodal diagnostic quantities 680 // --------------------------------------------------------------------------- 681 // Geometric factors 682 // -- Coordinate basis 683 CeedBasis basis_x; 684 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS_LOBATTO, 685 &basis_x); 686 // -- QFunction 687 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, 688 problem_data.setup_geo_loc, &qf_setup_geo); 689 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*dim, CEED_EVAL_GRAD); 690 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 691 CeedQFunctionAddOutput(qf_setup_geo, "geo_data", geo_data_size, CEED_EVAL_NONE); 692 // -- Operator 693 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, 694 CEED_QFUNCTION_NONE, &op_setup_geo); 695 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, 696 basis_x, CEED_VECTOR_ACTIVE); 697 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, 698 basis_x, CEED_VECTOR_NONE); 699 CeedOperatorSetField(op_setup_geo, "geo_data", 700 data[fine_level]->elem_restr_geo_data_diagnostic_i, 701 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 702 // -- Compute the quadrature data 703 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data_diagnostic, 704 CEED_REQUEST_IMMEDIATE); 705 // -- Cleanup 706 CeedBasisDestroy(&basis_x); 707 CeedQFunctionDestroy(&qf_setup_geo); 708 CeedOperatorDestroy(&op_setup_geo); 709 710 // Diagnostic quantities 711 // -- QFunction 712 CeedQFunctionCreateInterior(ceed, 1, problem_data.diagnostic, 713 problem_data.diagnostic_loc, &qf_diagnostic); 714 CeedQFunctionAddInput(qf_diagnostic, "u", num_comp_u, CEED_EVAL_INTERP); 715 CeedQFunctionAddInput(qf_diagnostic, "du", num_comp_u*dim, CEED_EVAL_GRAD); 716 CeedQFunctionAddInput(qf_diagnostic, "geo_data", geo_data_size, CEED_EVAL_NONE); 717 CeedQFunctionAddOutput(qf_diagnostic, "diagnostic", num_comp_u + num_comp_d, 718 CEED_EVAL_NONE); 719 CeedQFunctionSetContext(qf_diagnostic, phys_ctx); 720 // -- Operator 721 CeedOperatorCreate(ceed, qf_diagnostic, CEED_QFUNCTION_NONE, 722 CEED_QFUNCTION_NONE, &op_diagnostic); 723 CeedOperatorSetField(op_diagnostic, "u", data[fine_level]->elem_restr_u, 724 data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE); 725 CeedOperatorSetField(op_diagnostic, "du", data[fine_level]->elem_restr_u, 726 data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE); 727 CeedOperatorSetField(op_diagnostic, "geo_data", 728 data[fine_level]->elem_restr_geo_data_diagnostic_i, 729 CEED_BASIS_COLLOCATED, data[fine_level]->geo_data_diagnostic); 730 CeedOperatorSetField(op_diagnostic, "diagnostic", 731 data[fine_level]->elem_restr_diagnostic, 732 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 733 // -- Save libCEED data 734 data[fine_level]->qf_diagnostic = qf_diagnostic; 735 data[fine_level]->op_diagnostic = op_diagnostic; 736 737 // --------------------------------------------------------------------------- 738 // Cleanup 739 // --------------------------------------------------------------------------- 740 CeedVectorDestroy(&x_coord); 741 742 PetscFunctionReturn(0); 743 }; 744 745 // Set up libCEED multigrid level for a given degree 746 // Prolongation and Restriction are between level and level+1 747 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx, 748 ProblemData problem_data, PetscInt level, 749 PetscInt num_comp_u, PetscInt U_g_size, 750 PetscInt U_loc_size, CeedVector fine_mult, 751 CeedData *data) { 752 PetscErrorCode ierr; 753 CeedInt fine_level = app_ctx->num_levels - 1; 754 CeedInt P = app_ctx->level_degrees[level] + 1; 755 CeedInt Q = app_ctx->level_degrees[fine_level] + 1 + app_ctx->q_extra; 756 CeedInt dim; 757 CeedOperator op_jacobian, op_prolong, op_restrict; 758 759 PetscFunctionBeginUser; 760 761 ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr); 762 763 // --------------------------------------------------------------------------- 764 // libCEED restrictions 765 // --------------------------------------------------------------------------- 766 // -- Solution restriction 767 ierr = CreateRestrictionFromPlex(ceed, dm, P, 0, 0, 0, 768 &data[level]->elem_restr_u); 769 CHKERRQ(ierr); 770 771 // --------------------------------------------------------------------------- 772 // libCEED bases 773 // --------------------------------------------------------------------------- 774 // -- Solution basis 775 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, 776 problem_data.quadrature_mode, 777 &data[level]->basis_u); 778 779 // --------------------------------------------------------------------------- 780 // Persistent libCEED vectors 781 // --------------------------------------------------------------------------- 782 CeedVectorCreate(ceed, U_loc_size, &data[level]->x_ceed); 783 CeedVectorCreate(ceed, U_loc_size, &data[level]->y_ceed); 784 785 // --------------------------------------------------------------------------- 786 // Coarse Grid, Prolongation, and Restriction Operators 787 // --------------------------------------------------------------------------- 788 // Create the Operators that compute the prolongation and 789 // restriction between the p-multigrid levels and the coarse grid eval. 790 // --------------------------------------------------------------------------- 791 CeedOperatorMultigridLevelCreate(data[level+1]->op_jacobian, fine_mult, 792 data[level]->elem_restr_u, data[level]->basis_u, 793 &op_jacobian, &op_prolong, &op_restrict); 794 795 // -- Save libCEED data 796 data[level]->op_jacobian = op_jacobian; 797 data[level+1]->op_prolong = op_prolong; 798 data[level+1]->op_restrict = op_restrict; 799 800 PetscFunctionReturn(0); 801 }; 802