1 // Copyright (c) 2017-2025, 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 #include "../include/libceedsetup.h" 9 10 #include <stdio.h> 11 12 #include "../include/libceedsetup.h" 13 #include "../include/petscutils.h" 14 15 // ----------------------------------------------------------------------------- 16 // Destroy libCEED operator objects 17 // ----------------------------------------------------------------------------- 18 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) { 19 PetscFunctionBeginUser; 20 CeedVectorDestroy(&data->q_data); 21 CeedVectorDestroy(&data->x_ceed); 22 CeedVectorDestroy(&data->y_ceed); 23 CeedBasisDestroy(&data->basis_x); 24 CeedBasisDestroy(&data->basis_u); 25 CeedElemRestrictionDestroy(&data->elem_restr_u); 26 CeedElemRestrictionDestroy(&data->elem_restr_x); 27 CeedElemRestrictionDestroy(&data->elem_restr_u_i); 28 CeedElemRestrictionDestroy(&data->elem_restr_qd_i); 29 CeedQFunctionDestroy(&data->qf_apply); 30 CeedOperatorDestroy(&data->op_apply); 31 if (i > 0) { 32 CeedOperatorDestroy(&data->op_prolong); 33 CeedOperatorDestroy(&data->op_restrict); 34 } 35 PetscCall(PetscFree(data)); 36 PetscFunctionReturn(PETSC_SUCCESS); 37 }; 38 39 // ----------------------------------------------------------------------------- 40 // Set up libCEED for a given degree 41 // ----------------------------------------------------------------------------- 42 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, CeedInt topo_dim, CeedInt q_extra, PetscInt num_comp_x, PetscInt num_comp_u, 43 PetscInt g_size, PetscInt xl_size, BPData bp_data, CeedData data, PetscBool setup_rhs, PetscBool is_fine_level, 44 CeedVector rhs_ceed, CeedVector *target) { 45 DM dm_coord; 46 Vec coords; 47 const PetscScalar *coord_array; 48 CeedBasis basis_x, basis_u; 49 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 50 CeedQFunction qf_setup_geo = NULL, qf_apply = NULL; 51 CeedOperator op_setup_geo, op_apply; 52 CeedVector x_coord, q_data, x_ceed, y_ceed; 53 PetscInt c_start, c_end, num_elem; 54 CeedInt num_qpts, q_data_size = bp_data.q_data_size; 55 CeedScalar R = 1; // radius of the sphere 56 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube 57 58 PetscFunctionBeginUser; 59 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 60 61 // CEED bases 62 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x)); 63 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 64 65 // CEED restrictions 66 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x)); 67 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u)); 68 69 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); 70 num_elem = c_end - c_start; 71 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 72 73 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i); 74 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i); 75 76 // Element coordinates 77 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 78 PetscCall(VecGetArrayRead(coords, &coord_array)); 79 80 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 81 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array); 82 PetscCall(VecRestoreArrayRead(coords, &coord_array)); 83 84 // Create the persistent vectors that will be needed in setup and apply 85 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data); 86 CeedVectorCreate(ceed, xl_size, &x_ceed); 87 CeedVectorCreate(ceed, xl_size, &y_ceed); 88 89 if (is_fine_level) { 90 // Create the QFunction that builds the context data 91 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo); 92 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 93 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD); 94 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 95 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 96 97 // Create the operator that builds the quadrature data 98 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 99 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 100 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 101 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 102 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 103 104 // Setup q_data 105 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 106 107 // Set up PDE operator 108 PetscBool is_interp = bp_data.in_mode == CEED_EVAL_INTERP; 109 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1; 110 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1; 111 112 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply); 113 if (bp_data.in_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) { 114 CeedQFunctionAddInput(qf_apply, "u", num_comp_u, CEED_EVAL_INTERP); 115 CeedQFunctionAddInput(qf_apply, "du", num_comp_u * topo_dim, CEED_EVAL_GRAD); 116 } else { 117 CeedQFunctionAddInput(qf_apply, is_interp ? "u" : "du", num_comp_u * in_scale, bp_data.in_mode); 118 } 119 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 120 if (bp_data.out_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) { 121 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u, CEED_EVAL_INTERP); 122 CeedQFunctionAddOutput(qf_apply, "dv", num_comp_u * topo_dim, CEED_EVAL_GRAD); 123 } else { 124 CeedQFunctionAddOutput(qf_apply, is_interp ? "v" : "dv", num_comp_u * out_scale, bp_data.out_mode); 125 } 126 127 // Create the mass or diff operator 128 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 129 if (bp_data.in_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) { 130 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 131 CeedOperatorSetField(op_apply, "du", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 132 } else { 133 CeedOperatorSetField(op_apply, is_interp ? "u" : "du", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 134 } 135 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 136 if (bp_data.out_mode == CEED_EVAL_INTERP + CEED_EVAL_GRAD) { 137 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 138 CeedOperatorSetField(op_apply, "dv", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 139 } else { 140 CeedOperatorSetField(op_apply, is_interp ? "v" : "dv", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 141 } 142 143 // Cleanup 144 CeedQFunctionDestroy(&qf_setup_geo); 145 CeedOperatorDestroy(&op_setup_geo); 146 } 147 148 // Set up RHS if needed 149 if (setup_rhs) { 150 CeedQFunction qf_setup_rhs; 151 CeedOperator op_setup_rhs; 152 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, target); 153 // Create the q-function that sets up the RHS and true solution 154 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 155 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 156 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 157 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE); 158 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 159 160 // Create the operator that builds the RHS and true solution 161 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 162 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 163 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 164 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, *target); 165 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 166 167 // Set up the libCEED context 168 CeedQFunctionContext ctx_rhs_setup; 169 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 170 CeedScalar rhs_setup_data[2] = {R, l}; 171 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 172 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 173 CeedQFunctionContextDestroy(&ctx_rhs_setup); 174 175 // Setup RHS and target 176 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 177 178 // Cleanup 179 CeedQFunctionDestroy(&qf_setup_rhs); 180 CeedOperatorDestroy(&op_setup_rhs); 181 } 182 // Cleanup 183 CeedVectorDestroy(&x_coord); 184 185 // Save libCEED data required for level 186 data->basis_x = basis_x; 187 data->basis_u = basis_u; 188 data->elem_restr_x = elem_restr_x; 189 data->elem_restr_u = elem_restr_u; 190 data->elem_restr_u_i = elem_restr_u_i; 191 data->elem_restr_qd_i = elem_restr_qd_i; 192 data->qf_apply = qf_apply; 193 data->op_apply = op_apply; 194 data->q_data = q_data; 195 data->x_ceed = x_ceed; 196 data->y_ceed = y_ceed; 197 data->q_data_size = q_data_size; 198 PetscFunctionReturn(PETSC_SUCCESS); 199 }; 200 201 // ----------------------------------------------------------------------------- 202 // Setup libCEED level transfer operator objects 203 // ----------------------------------------------------------------------------- 204 PetscErrorCode CeedLevelTransferSetup(DM dm, Ceed ceed, CeedInt level, CeedInt num_comp_u, CeedData *data, BPData bp_data, Vec fine_mult) { 205 PetscFunctionBeginUser; 206 // Restriction - Fine to corse 207 CeedOperator op_restrict; 208 // Interpolation - Corse to fine 209 CeedOperator op_prolong; 210 // Coarse grid operator 211 CeedOperator op_apply; 212 // Basis 213 CeedBasis basis_u; 214 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 215 216 // --------------------------------------------------------------------------- 217 // Coarse Grid, Prolongation, and Restriction Operators 218 // --------------------------------------------------------------------------- 219 // Create the Operators that compute the prolongation and 220 // restriction between the p-multigrid levels and the coarse grid eval. 221 // --------------------------------------------------------------------------- 222 // Place in libCEED array 223 PetscMemType m_mem_type; 224 PetscCall(VecReadP2C(fine_mult, &m_mem_type, data[level]->x_ceed)); 225 226 CeedOperatorMultigridLevelCreate(data[level]->op_apply, data[level]->x_ceed, data[level - 1]->elem_restr_u, basis_u, &op_apply, &op_prolong, 227 &op_restrict); 228 229 // Restore PETSc vector 230 PetscCall(VecReadC2P(data[level]->x_ceed, m_mem_type, fine_mult)); 231 PetscCall(VecZeroEntries(fine_mult)); 232 // -- Save libCEED data 233 data[level - 1]->op_apply = op_apply; 234 data[level]->op_prolong = op_prolong; 235 data[level]->op_restrict = op_restrict; 236 237 CeedBasisDestroy(&basis_u); 238 PetscFunctionReturn(PETSC_SUCCESS); 239 }; 240 241 PetscErrorCode SetupErrorOperator(DM dm, Ceed ceed, BPData bp_data, CeedInt topo_dim, PetscInt num_comp_x, PetscInt num_comp_u, 242 CeedOperator *op_error) { 243 DM dm_coord; 244 Vec coords; 245 const PetscScalar *coord_array; 246 CeedBasis basis_x, basis_u; 247 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 248 CeedQFunction qf_setup_geo, qf_setup_rhs, qf_error; 249 CeedOperator op_setup_geo, op_setup_rhs; 250 CeedVector x_coord, q_data, target, rhs; 251 PetscInt c_start, c_end, num_elem; 252 CeedInt num_qpts, q_data_size = bp_data.q_data_size; 253 CeedScalar R = 1; // radius of the sphere 254 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube 255 256 PetscFunctionBeginUser; 257 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 258 259 // CEED bases 260 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x)); 261 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 262 263 // CEED restrictions 264 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x)); 265 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u)); 266 267 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); 268 num_elem = c_end - c_start; 269 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 270 271 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i); 272 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i); 273 274 // Element coordinates 275 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 276 PetscCall(VecGetArrayRead(coords, &coord_array)); 277 278 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 279 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array); 280 PetscCall(VecRestoreArrayRead(coords, &coord_array)); 281 282 // Create the persistent vectors that will be needed in setup and apply 283 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data); 284 285 // Create the QFunction that builds the context data 286 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo); 287 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 288 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD); 289 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 290 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 291 292 // Create the operator that builds the quadrature data 293 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 294 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 295 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 296 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 297 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 298 299 // Setup q_data 300 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 301 302 // Set up target vector 303 CeedElemRestrictionCreateVector(elem_restr_u, &rhs, NULL); 304 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, &target); 305 // Create the q-function that sets up the RHS and true solution 306 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 307 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 308 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 309 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE); 310 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 311 312 // Create the operator that builds the RHS and true solution 313 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 314 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 315 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 316 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, target); 317 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 318 319 // Set up the libCEED context 320 CeedQFunctionContext ctx_rhs_setup; 321 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 322 CeedScalar rhs_setup_data[2] = {R, l}; 323 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 324 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 325 CeedQFunctionContextDestroy(&ctx_rhs_setup); 326 327 // Setup RHS and target 328 CeedOperatorApply(op_setup_rhs, x_coord, rhs, CEED_REQUEST_IMMEDIATE); 329 330 // Set up error operator 331 // Create the error QFunction 332 CeedQFunctionCreateInterior(ceed, 1, bp_data.error, bp_data.error_loc, &qf_error); 333 CeedQFunctionAddInput(qf_error, "u", num_comp_u, CEED_EVAL_INTERP); 334 CeedQFunctionAddInput(qf_error, "true_soln", num_comp_u, CEED_EVAL_NONE); 335 CeedQFunctionAddInput(qf_error, "qdata", q_data_size, CEED_EVAL_NONE); 336 CeedQFunctionAddOutput(qf_error, "error", num_comp_u, CEED_EVAL_INTERP); 337 338 // Create the error operator 339 CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_error); 340 CeedOperatorSetField(*op_error, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 341 CeedOperatorSetField(*op_error, "true_soln", elem_restr_u_i, CEED_BASIS_NONE, target); 342 CeedOperatorSetField(*op_error, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 343 CeedOperatorSetField(*op_error, "error", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 344 345 // Cleanup 346 CeedQFunctionDestroy(&qf_setup_rhs); 347 CeedOperatorDestroy(&op_setup_rhs); 348 CeedQFunctionDestroy(&qf_setup_geo); 349 CeedOperatorDestroy(&op_setup_geo); 350 CeedQFunctionDestroy(&qf_error); 351 CeedVectorDestroy(&x_coord); 352 CeedVectorDestroy(&rhs); 353 CeedVectorDestroy(&target); 354 CeedVectorDestroy(&q_data); 355 CeedElemRestrictionDestroy(&elem_restr_u_i); 356 CeedElemRestrictionDestroy(&elem_restr_qd_i); 357 CeedElemRestrictionDestroy(&elem_restr_x); 358 CeedElemRestrictionDestroy(&elem_restr_u); 359 CeedBasisDestroy(&basis_x); 360 CeedBasisDestroy(&basis_u); 361 362 PetscFunctionReturn(PETSC_SUCCESS); 363 } 364