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