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 geo_data returns dXdx_i,j and w * det. 217 // --------------------------------------------------------------------------- 218 // -- QFunction 219 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, problem_data.setup_geo_loc, &qf_setup_geo); 220 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * dim, CEED_EVAL_GRAD); 221 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 222 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 223 // -- Operator 224 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_geo); 225 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, data[fine_level]->basis_x, CEED_VECTOR_ACTIVE); 226 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, data[fine_level]->basis_x, CEED_VECTOR_NONE); 227 CeedOperatorSetField(op_setup_geo, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 228 // -- Compute the quadrature data 229 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data, CEED_REQUEST_IMMEDIATE); 230 // -- Cleanup 231 CeedQFunctionDestroy(&qf_setup_geo); 232 CeedOperatorDestroy(&op_setup_geo); 233 234 // --------------------------------------------------------------------------- 235 // Local residual evaluator 236 // --------------------------------------------------------------------------- 237 // Create the QFunction and Operator that computes the residual of the non-linear PDE. 238 // --------------------------------------------------------------------------- 239 // -- QFunction 240 CeedQFunctionCreateInterior(ceed, 1, problem_data.residual, problem_data.residual_loc, &qf_residual); 241 CeedQFunctionAddInput(qf_residual, "du", num_comp_u * dim, CEED_EVAL_GRAD); 242 CeedQFunctionAddInput(qf_residual, "qdata", q_data_size, CEED_EVAL_NONE); 243 CeedQFunctionAddOutput(qf_residual, "dv", num_comp_u * dim, CEED_EVAL_GRAD); 244 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 245 CeedQFunctionAddOutput(qf_residual, problem_data.field_names[i], problem_data.field_sizes[i], CEED_EVAL_NONE); 246 } 247 CeedQFunctionSetContext(qf_residual, phys_ctx); 248 // -- Operator 249 CeedOperatorCreate(ceed, qf_residual, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_residual); 250 CeedOperatorSetField(op_residual, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 251 CeedOperatorSetField(op_residual, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 252 CeedOperatorSetField(op_residual, "dv", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 253 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 254 CeedOperatorSetField(op_residual, problem_data.field_names[i], data[fine_level]->elem_restr_stored_fields_i[i], CEED_BASIS_COLLOCATED, 255 data[fine_level]->stored_fields[i]); 256 } 257 // -- Save libCEED data 258 data[fine_level]->qf_residual = qf_residual; 259 data[fine_level]->op_residual = op_residual; 260 261 // --------------------------------------------------------------------------- 262 // Jacobian evaluator 263 // --------------------------------------------------------------------------- 264 // Create the QFunction and Operator that computes the action of the Jacobian for each linear solve. 265 // --------------------------------------------------------------------------- 266 // -- QFunction 267 CeedQFunctionCreateInterior(ceed, 1, problem_data.jacobian, problem_data.jacobian_loc, &qf_jacobian); 268 CeedQFunctionAddInput(qf_jacobian, "delta du", num_comp_u * dim, CEED_EVAL_GRAD); 269 CeedQFunctionAddInput(qf_jacobian, "qdata", q_data_size, CEED_EVAL_NONE); 270 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 271 CeedQFunctionAddInput(qf_jacobian, problem_data.field_names[i], problem_data.field_sizes[i], CEED_EVAL_NONE); 272 } 273 CeedQFunctionAddOutput(qf_jacobian, "delta dv", num_comp_u * dim, CEED_EVAL_GRAD); 274 CeedQFunctionSetContext(qf_jacobian, phys_ctx); 275 // -- Operator 276 CeedOperatorCreate(ceed, qf_jacobian, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_jacobian); 277 CeedOperatorSetField(op_jacobian, "delta du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 278 CeedOperatorSetField(op_jacobian, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 279 CeedOperatorSetField(op_jacobian, "delta dv", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 280 for (CeedInt i = 0; i < problem_data.number_fields_stored; i++) { 281 CeedOperatorSetField(op_jacobian, problem_data.field_names[i], data[fine_level]->elem_restr_stored_fields_i[i], CEED_BASIS_COLLOCATED, 282 data[fine_level]->stored_fields[i]); 283 } 284 // -- Save libCEED data 285 data[fine_level]->qf_jacobian = qf_jacobian; 286 data[fine_level]->op_jacobian = op_jacobian; 287 288 // --------------------------------------------------------------------------- 289 // Traction boundary conditions, if needed 290 // --------------------------------------------------------------------------- 291 if (app_ctx->bc_traction_count > 0) { 292 // -- Setup 293 DMLabel domain_label; 294 PetscCall(DMGetLabel(dm, "Face Sets", &domain_label)); 295 PetscCall(DMGetDimension(dm, &dim)); 296 297 // -- Basis 298 CeedInt height = 1; 299 CeedBasis basis_x_face, basis_u_face; 300 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_x, 2, Q, problem_data.quadrature_mode, &basis_x_face); 301 CeedBasisCreateTensorH1Lagrange(ceed, dim - height, num_comp_u, P, Q, problem_data.quadrature_mode, &basis_u_face); 302 // -- QFunction 303 CeedQFunction qf_traction; 304 CeedQFunctionContext traction_ctx; 305 CeedQFunctionCreateInterior(ceed, 1, SetupTractionBCs, SetupTractionBCs_loc, &qf_traction); 306 CeedQFunctionContextCreate(ceed, &traction_ctx); 307 CeedQFunctionSetContext(qf_traction, traction_ctx); 308 CeedQFunctionAddInput(qf_traction, "dx", num_comp_x * (num_comp_x - height), CEED_EVAL_GRAD); 309 CeedQFunctionAddInput(qf_traction, "weight", 1, CEED_EVAL_WEIGHT); 310 CeedQFunctionAddOutput(qf_traction, "v", num_comp_u, CEED_EVAL_INTERP); 311 312 // -- Compute contribution on each boundary face 313 for (CeedInt i = 0; i < app_ctx->bc_traction_count; i++) { 314 CeedElemRestriction elem_restr_x_face, elem_restr_u_face; 315 CeedOperator op_traction; 316 CeedQFunctionContextSetData(traction_ctx, CEED_MEM_HOST, CEED_USE_POINTER, 3 * sizeof(CeedScalar), app_ctx->bc_traction_vector[i]); 317 // Setup restriction 318 PetscCall( 319 GetRestrictionForDomain(ceed, dm, 1, domain_label, app_ctx->bc_traction_faces[i], Q, 0, &elem_restr_u_face, &elem_restr_x_face, NULL)); 320 // ---- Create boundary Operator 321 CeedOperatorCreate(ceed, qf_traction, NULL, NULL, &op_traction); 322 CeedOperatorSetField(op_traction, "dx", elem_restr_x_face, basis_x_face, CEED_VECTOR_ACTIVE); 323 CeedOperatorSetField(op_traction, "weight", CEED_ELEMRESTRICTION_NONE, basis_x_face, CEED_VECTOR_NONE); 324 CeedOperatorSetField(op_traction, "v", elem_restr_u_face, basis_u_face, CEED_VECTOR_ACTIVE); 325 // ---- Compute traction on face 326 CeedOperatorApplyAdd(op_traction, x_coord, neumann_ceed, CEED_REQUEST_IMMEDIATE); 327 // ---- Cleanup 328 CeedElemRestrictionDestroy(&elem_restr_x_face); 329 CeedElemRestrictionDestroy(&elem_restr_u_face); 330 CeedOperatorDestroy(&op_traction); 331 } 332 // -- Cleanup 333 CeedBasisDestroy(&basis_x_face); 334 CeedBasisDestroy(&basis_u_face); 335 CeedQFunctionDestroy(&qf_traction); 336 CeedQFunctionContextDestroy(&traction_ctx); 337 } 338 339 // --------------------------------------------------------------------------- 340 // Forcing term, if needed 341 // --------------------------------------------------------------------------- 342 // Create the QFunction and Operator that computes the forcing term (RHS) for the non-linear PDE. 343 // --------------------------------------------------------------------------- 344 if (forcing_choice != FORCE_NONE) { 345 CeedQFunction qf_setup_force; 346 CeedOperator op_setup_force; 347 348 // -- QFunction 349 CeedQFunctionCreateInterior(ceed, 1, forcing_options[forcing_choice].setup_forcing, forcing_options[forcing_choice].setup_forcing_loc, 350 &qf_setup_force); 351 CeedQFunctionAddInput(qf_setup_force, "x", num_comp_x, CEED_EVAL_INTERP); 352 CeedQFunctionAddInput(qf_setup_force, "qdata", q_data_size, CEED_EVAL_NONE); 353 CeedQFunctionAddOutput(qf_setup_force, "force", num_comp_u, CEED_EVAL_INTERP); 354 if (forcing_choice == FORCE_MMS) { 355 CeedQFunctionSetContext(qf_setup_force, phys_ctx); 356 } else { 357 CeedQFunctionContext ctxForcing; 358 CeedQFunctionContextCreate(ceed, &ctxForcing); 359 CeedQFunctionContextSetData(ctxForcing, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*app_ctx->forcing_vector), app_ctx->forcing_vector); 360 CeedQFunctionSetContext(qf_setup_force, ctxForcing); 361 CeedQFunctionContextDestroy(&ctxForcing); 362 } 363 // -- Operator 364 CeedOperatorCreate(ceed, qf_setup_force, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_force); 365 CeedOperatorSetField(op_setup_force, "x", data[fine_level]->elem_restr_x, data[fine_level]->basis_x, CEED_VECTOR_ACTIVE); 366 CeedOperatorSetField(op_setup_force, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 367 CeedOperatorSetField(op_setup_force, "force", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 368 // -- Compute forcing term 369 CeedOperatorApply(op_setup_force, x_coord, force_ceed, CEED_REQUEST_IMMEDIATE); 370 // -- Cleanup 371 CeedQFunctionDestroy(&qf_setup_force); 372 CeedOperatorDestroy(&op_setup_force); 373 } 374 375 // --------------------------------------------------------------------------- 376 // True solution, for MMS 377 // --------------------------------------------------------------------------- 378 // Create the QFunction and Operator that computes the true solution at the mesh nodes for validation with the manufactured solution. 379 // --------------------------------------------------------------------------- 380 if (problem_data.true_soln) { 381 CeedScalar *true_array; 382 const CeedScalar *mult_array; 383 CeedVector mult_vec; 384 CeedBasis basis_x_true; 385 CeedQFunction qf_true; 386 CeedOperator op_true; 387 388 // -- Solution vector 389 CeedVectorCreate(ceed, U_loc_size, &(data[fine_level]->true_soln)); 390 // -- Basis 391 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, P, CEED_GAUSS_LOBATTO, &basis_x_true); 392 // QFunction 393 CeedQFunctionCreateInterior(ceed, 1, problem_data.true_soln, problem_data.true_soln_loc, &qf_true); 394 CeedQFunctionAddInput(qf_true, "x", num_comp_x, CEED_EVAL_INTERP); 395 CeedQFunctionAddOutput(qf_true, "true solution", num_comp_u, CEED_EVAL_NONE); 396 // Operator 397 CeedOperatorCreate(ceed, qf_true, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_true); 398 CeedOperatorSetField(op_true, "x", data[fine_level]->elem_restr_x, basis_x_true, CEED_VECTOR_ACTIVE); 399 CeedOperatorSetField(op_true, "true solution", data[fine_level]->elem_restr_u, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 400 // -- Compute true solution 401 CeedOperatorApply(op_true, x_coord, data[fine_level]->true_soln, CEED_REQUEST_IMMEDIATE); 402 // -- Multiplicity calculation 403 CeedElemRestrictionCreateVector(data[fine_level]->elem_restr_u, &mult_vec, NULL); 404 CeedVectorSetValue(mult_vec, 0.); 405 CeedElemRestrictionGetMultiplicity(data[fine_level]->elem_restr_u, mult_vec); 406 // -- Multiplicity correction 407 CeedVectorGetArray(data[fine_level]->true_soln, CEED_MEM_HOST, &true_array); 408 CeedVectorGetArrayRead(mult_vec, CEED_MEM_HOST, &mult_array); 409 for (CeedInt i = 0; i < U_loc_size; i++) true_array[i] /= mult_array[i]; 410 CeedVectorRestoreArray(data[fine_level]->true_soln, &true_array); 411 CeedVectorRestoreArrayRead(mult_vec, &mult_array); 412 // -- Cleanup 413 CeedVectorDestroy(&mult_vec); 414 CeedBasisDestroy(&basis_x_true); 415 CeedQFunctionDestroy(&qf_true); 416 CeedOperatorDestroy(&op_true); 417 } 418 419 // --------------------------------------------------------------------------- 420 // Local energy computation 421 // --------------------------------------------------------------------------- 422 // Create the QFunction and Operator that computes the strain energy 423 // --------------------------------------------------------------------------- 424 // -- QFunction 425 CeedQFunctionCreateInterior(ceed, 1, problem_data.energy, problem_data.energy_loc, &qf_energy); 426 CeedQFunctionAddInput(qf_energy, "du", num_comp_u * dim, CEED_EVAL_GRAD); 427 CeedQFunctionAddInput(qf_energy, "qdata", q_data_size, CEED_EVAL_NONE); 428 CeedQFunctionAddOutput(qf_energy, "energy", num_comp_e, CEED_EVAL_INTERP); 429 CeedQFunctionSetContext(qf_energy, phys_ctx); 430 // -- Operator 431 CeedOperatorCreate(ceed, qf_energy, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_energy); 432 CeedOperatorSetField(op_energy, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_u, CEED_VECTOR_ACTIVE); 433 CeedOperatorSetField(op_energy, "qdata", data[fine_level]->elem_restr_geo_data_i, CEED_BASIS_COLLOCATED, data[fine_level]->geo_data); 434 CeedOperatorSetField(op_energy, "energy", data[fine_level]->elem_restr_energy, data[fine_level]->basis_energy, CEED_VECTOR_ACTIVE); 435 // -- Save libCEED data 436 data[fine_level]->qf_energy = qf_energy; 437 data[fine_level]->op_energy = op_energy; 438 439 // --------------------------------------------------------------------------- 440 // Diagnostic value computation 441 // --------------------------------------------------------------------------- 442 // Create the QFunction and Operator that computes nodal diagnostic quantities 443 // --------------------------------------------------------------------------- 444 // Geometric factors 445 // -- Coordinate basis 446 CeedBasis basis_x; 447 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q, CEED_GAUSS_LOBATTO, &basis_x); 448 // -- QFunction 449 CeedQFunctionCreateInterior(ceed, 1, problem_data.setup_geo, problem_data.setup_geo_loc, &qf_setup_geo); 450 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * dim, CEED_EVAL_GRAD); 451 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 452 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 453 // -- Operator 454 CeedOperatorCreate(ceed, qf_setup_geo, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup_geo); 455 CeedOperatorSetField(op_setup_geo, "dx", data[fine_level]->elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 456 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 457 CeedOperatorSetField(op_setup_geo, "qdata", data[fine_level]->elem_restr_geo_data_diagnostic_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 458 // -- Compute the quadrature data 459 CeedOperatorApply(op_setup_geo, x_coord, data[fine_level]->geo_data_diagnostic, CEED_REQUEST_IMMEDIATE); 460 // -- Cleanup 461 CeedBasisDestroy(&basis_x); 462 CeedQFunctionDestroy(&qf_setup_geo); 463 CeedOperatorDestroy(&op_setup_geo); 464 465 // Diagnostic quantities 466 // -- QFunction 467 CeedQFunctionCreateInterior(ceed, 1, problem_data.diagnostic, problem_data.diagnostic_loc, &qf_diagnostic); 468 CeedQFunctionAddInput(qf_diagnostic, "u", num_comp_u, CEED_EVAL_INTERP); 469 CeedQFunctionAddInput(qf_diagnostic, "du", num_comp_u * dim, CEED_EVAL_GRAD); 470 CeedQFunctionAddInput(qf_diagnostic, "qdata", q_data_size, CEED_EVAL_NONE); 471 CeedQFunctionAddOutput(qf_diagnostic, "diagnostic values", num_comp_u + num_comp_d, CEED_EVAL_NONE); 472 CeedQFunctionSetContext(qf_diagnostic, phys_ctx); 473 // -- Operator 474 CeedOperatorCreate(ceed, qf_diagnostic, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_diagnostic); 475 CeedOperatorSetField(op_diagnostic, "u", data[fine_level]->elem_restr_u, data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE); 476 CeedOperatorSetField(op_diagnostic, "du", data[fine_level]->elem_restr_u, data[fine_level]->basis_diagnostic, CEED_VECTOR_ACTIVE); 477 CeedOperatorSetField(op_diagnostic, "qdata", data[fine_level]->elem_restr_geo_data_diagnostic_i, CEED_BASIS_COLLOCATED, 478 data[fine_level]->geo_data_diagnostic); 479 CeedOperatorSetField(op_diagnostic, "diagnostic values", data[fine_level]->elem_restr_diagnostic, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 480 // -- Save libCEED data 481 data[fine_level]->qf_diagnostic = qf_diagnostic; 482 data[fine_level]->op_diagnostic = op_diagnostic; 483 484 // --------------------------------------------------------------------------- 485 // Cleanup 486 // --------------------------------------------------------------------------- 487 CeedVectorDestroy(&x_coord); 488 489 PetscFunctionReturn(0); 490 }; 491 492 // Set up libCEED multigrid level for a given degree 493 // Prolongation and Restriction are between level and level+1 494 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx, ProblemData problem_data, PetscInt level, PetscInt num_comp_u, PetscInt U_g_size, 495 PetscInt U_loc_size, CeedVector fine_mult, CeedData *data) { 496 CeedInt fine_level = app_ctx->num_levels - 1; 497 CeedInt P = app_ctx->level_degrees[level] + 1; 498 CeedInt Q = app_ctx->level_degrees[fine_level] + 1 + app_ctx->q_extra; 499 CeedInt dim; 500 CeedOperator op_jacobian, op_prolong, op_restrict; 501 502 PetscFunctionBeginUser; 503 504 PetscCall(DMGetDimension(dm, &dim)); 505 506 // --------------------------------------------------------------------------- 507 // libCEED restrictions 508 // --------------------------------------------------------------------------- 509 // -- Solution restriction 510 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &data[level]->elem_restr_u)); 511 512 // --------------------------------------------------------------------------- 513 // libCEED bases 514 // --------------------------------------------------------------------------- 515 // -- Solution basis 516 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_u, P, Q, problem_data.quadrature_mode, &data[level]->basis_u); 517 518 // --------------------------------------------------------------------------- 519 // Persistent libCEED vectors 520 // --------------------------------------------------------------------------- 521 CeedVectorCreate(ceed, U_loc_size, &data[level]->x_ceed); 522 CeedVectorCreate(ceed, U_loc_size, &data[level]->y_ceed); 523 524 // --------------------------------------------------------------------------- 525 // Coarse Grid, Prolongation, and Restriction Operators 526 // --------------------------------------------------------------------------- 527 // Create the Operators that compute the prolongation and restriction between the p-multigrid levels and the coarse grid eval. 528 // --------------------------------------------------------------------------- 529 CeedOperatorMultigridLevelCreate(data[level + 1]->op_jacobian, fine_mult, data[level]->elem_restr_u, data[level]->basis_u, &op_jacobian, 530 &op_prolong, &op_restrict); 531 532 // -- Save libCEED data 533 data[level]->op_jacobian = op_jacobian; 534 data[level + 1]->op_prolong = op_prolong; 535 data[level + 1]->op_restrict = op_restrict; 536 537 PetscFunctionReturn(0); 538 }; 539