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