1 #include <stdio.h> 2 #include "../include/libceedsetup.h" 3 #include "../include/petscutils.h" 4 5 // ----------------------------------------------------------------------------- 6 // Destroy libCEED operator objects 7 // ----------------------------------------------------------------------------- 8 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) { 9 int ierr; 10 11 PetscFunctionBeginUser; 12 CeedVectorDestroy(&data->q_data); 13 CeedVectorDestroy(&data->x_ceed); 14 CeedVectorDestroy(&data->y_ceed); 15 CeedBasisDestroy(&data->basis_x); 16 CeedBasisDestroy(&data->basis_u); 17 CeedElemRestrictionDestroy(&data->elem_restr_u); 18 CeedElemRestrictionDestroy(&data->elem_restr_x); 19 CeedElemRestrictionDestroy(&data->elem_restr_u_i); 20 CeedElemRestrictionDestroy(&data->elem_restr_qd_i); 21 CeedQFunctionDestroy(&data->qf_apply); 22 CeedOperatorDestroy(&data->op_apply); 23 if (i > 0) { 24 CeedOperatorDestroy(&data->op_prolong); 25 CeedBasisDestroy(&data->basis_c_to_f); 26 CeedOperatorDestroy(&data->op_restrict); 27 } 28 ierr = PetscFree(data); CHKERRQ(ierr); 29 30 PetscFunctionReturn(0); 31 }; 32 33 // ----------------------------------------------------------------------------- 34 // Set up libCEED for a given degree 35 // ----------------------------------------------------------------------------- 36 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, 37 CeedInt topo_dim, CeedInt q_extra, 38 PetscInt num_comp_x, PetscInt num_comp_u, 39 PetscInt g_size, PetscInt xl_size, 40 BPData bp_data, CeedData data, 41 PetscBool setup_rhs, CeedVector rhs_ceed, 42 CeedVector *target) { 43 int ierr; 44 DM dm_coord; 45 Vec coords; 46 const PetscScalar *coord_array; 47 CeedBasis basis_x, basis_u; 48 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 49 CeedQFunction qf_setup_geo, qf_apply; 50 CeedOperator op_setup_geo, op_apply; 51 CeedVector x_coord, q_data, x_ceed, y_ceed; 52 CeedInt num_qpts, c_start, c_end, num_elem, 53 q_data_size = bp_data.q_data_size; 54 CeedScalar R = 1, // radius of the sphere 55 l = 1.0/PetscSqrtReal(3.0); // half edge of the inscribed cube 56 57 <<<<<<< HEAD 58 PetscFunctionBeginUser; 59 // CEED bases 60 P = degree + 1; 61 Q = P + q_extra; 62 CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_u, P, Q, 63 bp_data.q_mode, 64 &basis_u); 65 CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_x, 2, Q, 66 bp_data.q_mode, 67 &basis_x); 68 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 69 70 // CEED restrictions 71 ======= 72 >>>>>>> 158419b6 (example/petsc: added CreateBasisFromPlex and tested with tensor basis) 73 ierr = DMSetCoordinateDim(dm, topo_dim); CHKERRQ(ierr); 74 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 75 76 // CEED bases 77 ierr = CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, &basis_x); CHKERRQ(ierr); 78 ierr = CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, &basis_u); CHKERRQ(ierr); 79 80 // CEED restrictions 81 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 82 CHKERRQ(ierr); 83 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x); 84 CHKERRQ(ierr); 85 ierr = CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u); 86 CHKERRQ(ierr); 87 88 ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); CHKERRQ(ierr); 89 num_elem = c_end - c_start; 90 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 91 92 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, 93 num_comp_u*num_elem*num_qpts, 94 CEED_STRIDES_BACKEND, &elem_restr_u_i); 95 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, 96 q_data_size*num_elem*num_qpts, 97 CEED_STRIDES_BACKEND, &elem_restr_qd_i); 98 99 // Element coordinates 100 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 101 ierr = VecGetArrayRead(coords, &coord_array); CHKERRQ(ierr); 102 103 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 104 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 105 (PetscScalar *)coord_array); 106 ierr = VecRestoreArrayRead(coords, &coord_array); 107 108 // Create the persistent vectors that will be needed in setup and apply 109 CeedVectorCreate(ceed, q_data_size*num_elem*num_qpts, &q_data); 110 CeedVectorCreate(ceed, xl_size, &x_ceed); 111 CeedVectorCreate(ceed, xl_size, &y_ceed); 112 113 // Create the QFunction that builds the context data 114 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, 115 &qf_setup_geo); 116 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 117 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*topo_dim, CEED_EVAL_GRAD); 118 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 119 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 120 121 // Create the operator that builds the quadrature data 122 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 123 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, 124 CEED_VECTOR_ACTIVE); 125 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, 126 CEED_VECTOR_ACTIVE); 127 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, 128 CEED_VECTOR_NONE); 129 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, 130 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 131 132 // Setup q_data 133 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 134 135 // Set up PDE operator 136 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1; 137 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1; 138 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, 139 &qf_apply); 140 CeedQFunctionAddInput(qf_apply, "u", num_comp_u*in_scale, bp_data.in_mode); 141 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 142 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u*out_scale, bp_data.out_mode); 143 144 // Create the mass or diff operator 145 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 146 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 147 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_COLLOCATED, 148 q_data); 149 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 150 151 // Set up RHS if needed 152 if (setup_rhs) { 153 CeedQFunction qf_setup_rhs; 154 CeedOperator op_setup_rhs; 155 CeedVectorCreate(ceed, num_elem*num_qpts*num_comp_u, target); 156 157 // Create the q-function that sets up the RHS and true solution 158 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, 159 &qf_setup_rhs); 160 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 161 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 162 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, 163 CEED_EVAL_NONE); 164 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 165 166 // Create the operator that builds the RHS and true solution 167 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 168 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, 169 CEED_VECTOR_ACTIVE); 170 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, 171 CEED_BASIS_COLLOCATED, q_data); 172 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, 173 CEED_BASIS_COLLOCATED, *target); 174 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, 175 CEED_VECTOR_ACTIVE); 176 177 // Set up the libCEED context 178 CeedQFunctionContext ctx_rhs_setup; 179 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 180 CeedScalar rhs_setup_data[2] = {R, l}; 181 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, 182 sizeof rhs_setup_data, &rhs_setup_data); 183 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 184 CeedQFunctionContextDestroy(&ctx_rhs_setup); 185 186 // Setup RHS and target 187 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 188 189 // Cleanup 190 CeedQFunctionDestroy(&qf_setup_rhs); 191 CeedOperatorDestroy(&op_setup_rhs); 192 } 193 194 // Cleanup 195 CeedQFunctionDestroy(&qf_setup_geo); 196 CeedOperatorDestroy(&op_setup_geo); 197 CeedVectorDestroy(&x_coord); 198 199 // Save libCEED data required for level 200 data->basis_x = basis_x; data->basis_u = basis_u; 201 data->elem_restr_x = elem_restr_x; 202 data->elem_restr_u = elem_restr_u; 203 data->elem_restr_u_i = elem_restr_u_i; 204 data->elem_restr_qd_i = elem_restr_qd_i; 205 data->qf_apply = qf_apply; 206 data->op_apply = op_apply; 207 data->q_data = q_data; 208 data->x_ceed = x_ceed; 209 data->y_ceed = y_ceed; 210 211 PetscFunctionReturn(0); 212 }; 213 214 // ----------------------------------------------------------------------------- 215 // Setup libCEED level transfer operator objects 216 // ----------------------------------------------------------------------------- 217 PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels, 218 CeedInt num_comp_u, CeedData *data, 219 CeedInt *level_degrees, 220 CeedQFunction qf_restrict, CeedQFunction qf_prolong) { 221 PetscFunctionBeginUser; 222 // Return early if num_levels=1 223 if (num_levels == 1) 224 PetscFunctionReturn(0); 225 226 // Set up each level 227 for (CeedInt i=1; i<num_levels; i++) { 228 // P coarse and P fine 229 CeedInt Pc = level_degrees[i-1] + 1; 230 CeedInt Pf = level_degrees[i] + 1; 231 232 // Restriction - Fine to corse 233 CeedBasis basis_c_to_f; 234 CeedOperator op_restrict; 235 236 // Basis 237 CeedBasisCreateTensorH1Lagrange(ceed, 3, num_comp_u, Pc, Pf, 238 CEED_GAUSS_LOBATTO, &basis_c_to_f); 239 240 // Create the restriction operator 241 CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, 242 CEED_QFUNCTION_NONE, &op_restrict); 243 CeedOperatorSetField(op_restrict, "input", data[i]->elem_restr_u, 244 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 245 CeedOperatorSetField(op_restrict, "output", data[i-1]->elem_restr_u, 246 basis_c_to_f, CEED_VECTOR_ACTIVE); 247 248 // Save libCEED data required for level 249 data[i]->basis_c_to_f = basis_c_to_f; 250 data[i]->op_restrict = op_restrict; 251 252 // Interpolation - Corse to fine 253 CeedOperator op_prolong; 254 255 // Create the prolongation operator 256 CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, 257 CEED_QFUNCTION_NONE, &op_prolong); 258 CeedOperatorSetField(op_prolong, "input", data[i-1]->elem_restr_u, 259 basis_c_to_f, CEED_VECTOR_ACTIVE); 260 CeedOperatorSetField(op_prolong, "output", data[i]->elem_restr_u, 261 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 262 263 // Save libCEED data required for level 264 data[i]->op_prolong = op_prolong; 265 } 266 267 PetscFunctionReturn(0); 268 }; 269 270 // ----------------------------------------------------------------------------- 271