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