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