1 #include "../include/libceedsetup.h" 2 3 #include <stdio.h> 4 5 #include "../include/petscutils.h" 6 7 // ----------------------------------------------------------------------------- 8 // Destroy libCEED operator objects 9 // ----------------------------------------------------------------------------- 10 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) { 11 PetscFunctionBeginUser; 12 13 CeedVectorDestroy(&data->q_data); 14 CeedVectorDestroy(&data->x_ceed); 15 CeedVectorDestroy(&data->y_ceed); 16 CeedBasisDestroy(&data->basis_x); 17 CeedBasisDestroy(&data->basis_u); 18 CeedElemRestrictionDestroy(&data->elem_restr_u); 19 CeedElemRestrictionDestroy(&data->elem_restr_x); 20 CeedElemRestrictionDestroy(&data->elem_restr_u_i); 21 CeedElemRestrictionDestroy(&data->elem_restr_qd_i); 22 CeedQFunctionDestroy(&data->qf_apply); 23 CeedOperatorDestroy(&data->op_apply); 24 if (i > 0) { 25 CeedOperatorDestroy(&data->op_prolong); 26 CeedOperatorDestroy(&data->op_restrict); 27 } 28 PetscCall(PetscFree(data)); 29 30 PetscFunctionReturn(PETSC_SUCCESS); 31 }; 32 33 // ----------------------------------------------------------------------------- 34 // Set up libCEED for a given degree 35 // ----------------------------------------------------------------------------- 36 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, CeedInt topo_dim, CeedInt q_extra, PetscInt num_comp_x, PetscInt num_comp_u, 37 PetscInt g_size, PetscInt xl_size, BPData bp_data, CeedData data, PetscBool setup_rhs, CeedVector rhs_ceed, 38 CeedVector *target) { 39 DM dm_coord; 40 Vec coords; 41 const PetscScalar *coord_array; 42 CeedBasis basis_x, basis_u; 43 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 44 CeedQFunction qf_setup_geo, qf_apply; 45 CeedOperator op_setup_geo, op_apply; 46 CeedVector x_coord, q_data, x_ceed, y_ceed; 47 CeedInt num_qpts, c_start, c_end, num_elem, q_data_size = bp_data.q_data_size; 48 CeedScalar R = 1; // radius of the sphere 49 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube 50 51 PetscFunctionBeginUser; 52 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 53 54 // CEED bases 55 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x)); 56 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 57 58 // CEED restrictions 59 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x)); 60 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u)); 61 62 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); 63 num_elem = c_end - c_start; 64 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 65 66 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i); 67 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i); 68 69 // Element coordinates 70 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 71 PetscCall(VecGetArrayRead(coords, &coord_array)); 72 73 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 74 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array); 75 PetscCall(VecRestoreArrayRead(coords, &coord_array)); 76 77 // Create the persistent vectors that will be needed in setup and apply 78 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data); 79 CeedVectorCreate(ceed, xl_size, &x_ceed); 80 CeedVectorCreate(ceed, xl_size, &y_ceed); 81 82 // Create the QFunction that builds the context data 83 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo); 84 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 85 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD); 86 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 87 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 88 89 // Create the operator that builds the quadrature data 90 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 91 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 92 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 93 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 94 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 95 96 // Setup q_data 97 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 98 99 // Set up PDE operator 100 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1; 101 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1; 102 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply); 103 CeedQFunctionAddInput(qf_apply, "u", num_comp_u * in_scale, bp_data.in_mode); 104 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 105 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u * out_scale, bp_data.out_mode); 106 107 // Create the mass or diff operator 108 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 109 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 110 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 111 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 112 113 // Set up RHS if needed 114 if (setup_rhs) { 115 CeedQFunction qf_setup_rhs; 116 CeedOperator op_setup_rhs; 117 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, target); 118 // Create the q-function that sets up the RHS and true solution 119 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 120 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 121 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 122 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE); 123 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 124 125 // Create the operator that builds the RHS and true solution 126 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 127 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 128 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 129 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, *target); 130 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 131 132 // Set up the libCEED context 133 CeedQFunctionContext ctx_rhs_setup; 134 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 135 CeedScalar rhs_setup_data[2] = {R, l}; 136 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 137 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 138 CeedQFunctionContextDestroy(&ctx_rhs_setup); 139 140 // Setup RHS and target 141 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 142 143 // Cleanup 144 CeedQFunctionDestroy(&qf_setup_rhs); 145 CeedOperatorDestroy(&op_setup_rhs); 146 } 147 148 // Cleanup 149 CeedQFunctionDestroy(&qf_setup_geo); 150 CeedOperatorDestroy(&op_setup_geo); 151 CeedVectorDestroy(&x_coord); 152 153 // Save libCEED data required for level 154 data->basis_x = basis_x; 155 data->basis_u = basis_u; 156 data->elem_restr_x = elem_restr_x; 157 data->elem_restr_u = elem_restr_u; 158 data->elem_restr_u_i = elem_restr_u_i; 159 data->elem_restr_qd_i = elem_restr_qd_i; 160 data->qf_apply = qf_apply; 161 data->op_apply = op_apply; 162 data->q_data = q_data; 163 data->x_ceed = x_ceed; 164 data->y_ceed = y_ceed; 165 data->q_data_size = q_data_size; 166 167 PetscFunctionReturn(PETSC_SUCCESS); 168 }; 169 170 // ----------------------------------------------------------------------------- 171 // Setup libCEED level transfer operator objects 172 // ----------------------------------------------------------------------------- 173 PetscErrorCode CeedLevelTransferSetup(DM dm, Ceed ceed, CeedInt level, CeedInt num_comp_u, CeedData *data, BPData bp_data, Vec fine_mult) { 174 PetscFunctionBeginUser; 175 // Restriction - Fine to corse 176 CeedOperator op_restrict; 177 // Interpolation - Corse to fine 178 CeedOperator op_prolong; 179 // Coarse grid operator 180 CeedOperator op_apply; 181 // Basis 182 CeedBasis basis_u; 183 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 184 185 // --------------------------------------------------------------------------- 186 // Coarse Grid, Prolongation, and Restriction Operators 187 // --------------------------------------------------------------------------- 188 // Create the Operators that compute the prolongation and 189 // restriction between the p-multigrid levels and the coarse grid eval. 190 // --------------------------------------------------------------------------- 191 // Place in libCEED array 192 const PetscScalar *m; 193 PetscMemType m_mem_type; 194 PetscCall(VecGetArrayReadAndMemType(fine_mult, &m, &m_mem_type)); 195 CeedVectorSetArray(data[level]->x_ceed, MemTypeP2C(m_mem_type), CEED_USE_POINTER, (CeedScalar *)m); 196 197 CeedOperatorMultigridLevelCreate(data[level]->op_apply, data[level]->x_ceed, data[level - 1]->elem_restr_u, basis_u, &op_apply, &op_prolong, 198 &op_restrict); 199 200 // Restore PETSc vector 201 CeedVectorTakeArray(data[level]->x_ceed, MemTypeP2C(m_mem_type), (CeedScalar **)&m); 202 PetscCall(VecRestoreArrayReadAndMemType(fine_mult, &m)); 203 PetscCall(VecZeroEntries(fine_mult)); 204 // -- Save libCEED data 205 data[level - 1]->op_apply = op_apply; 206 data[level]->op_prolong = op_prolong; 207 data[level]->op_restrict = op_restrict; 208 209 CeedBasisDestroy(&basis_u); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 }; 212 213 // ----------------------------------------------------------------------------- 214