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 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 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 <<<<<<< HEAD 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 ======= 75 //ierr = DMSetCoordinateDim(dm, topo_dim); CHKERRQ(ierr); 76 >>>>>>> 0fa86f50 (example/petsc: Added CreateDistributedDM in petscutils.c and some cleanup) 77 ierr = DMGetCoordinateDM(dm, &dm_coord); CHKERRQ(ierr); 78 79 // CEED bases 80 ierr = CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, &basis_x); CHKERRQ(ierr); 81 ierr = CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, &basis_u); CHKERRQ(ierr); 82 83 // CEED restrictions 84 ierr = DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL); 85 CHKERRQ(ierr); 86 ierr = CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x); 87 CHKERRQ(ierr); 88 ierr = CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u); 89 CHKERRQ(ierr); 90 91 ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); CHKERRQ(ierr); 92 num_elem = c_end - c_start; 93 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 94 95 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, 96 num_comp_u*num_elem*num_qpts, 97 CEED_STRIDES_BACKEND, &elem_restr_u_i); 98 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, 99 q_data_size*num_elem*num_qpts, 100 CEED_STRIDES_BACKEND, &elem_restr_qd_i); 101 102 // Element coordinates 103 ierr = DMGetCoordinatesLocal(dm, &coords); CHKERRQ(ierr); 104 ierr = VecGetArrayRead(coords, &coord_array); CHKERRQ(ierr); 105 106 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 107 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, 108 (PetscScalar *)coord_array); 109 ierr = VecRestoreArrayRead(coords, &coord_array); 110 111 // Create the persistent vectors that will be needed in setup and apply 112 CeedVectorCreate(ceed, q_data_size*num_elem*num_qpts, &q_data); 113 CeedVectorCreate(ceed, xl_size, &x_ceed); 114 CeedVectorCreate(ceed, xl_size, &y_ceed); 115 116 // Create the QFunction that builds the context data 117 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, 118 &qf_setup_geo); 119 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 120 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x*topo_dim, CEED_EVAL_GRAD); 121 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 122 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 123 124 // Create the operator that builds the quadrature data 125 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 126 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, 127 CEED_VECTOR_ACTIVE); 128 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, 129 CEED_VECTOR_ACTIVE); 130 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, 131 CEED_VECTOR_NONE); 132 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, 133 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 134 135 // Setup q_data 136 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 137 138 // Set up PDE operator 139 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1; 140 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1; 141 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, 142 &qf_apply); 143 CeedQFunctionAddInput(qf_apply, "u", num_comp_u*in_scale, bp_data.in_mode); 144 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 145 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u*out_scale, bp_data.out_mode); 146 147 // Create the mass or diff operator 148 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 149 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 150 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_COLLOCATED, 151 q_data); 152 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 153 154 // Set up RHS if needed 155 if (setup_rhs) { 156 CeedQFunction qf_setup_rhs; 157 CeedOperator op_setup_rhs; 158 CeedVectorCreate(ceed, num_elem*num_qpts*num_comp_u, target); 159 160 // Create the q-function that sets up the RHS and true solution 161 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, 162 &qf_setup_rhs); 163 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 164 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 165 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, 166 CEED_EVAL_NONE); 167 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 168 169 // Create the operator that builds the RHS and true solution 170 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 171 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, 172 CEED_VECTOR_ACTIVE); 173 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, 174 CEED_BASIS_COLLOCATED, q_data); 175 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, 176 CEED_BASIS_COLLOCATED, *target); 177 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, 178 CEED_VECTOR_ACTIVE); 179 180 // Set up the libCEED context 181 CeedQFunctionContext ctx_rhs_setup; 182 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 183 CeedScalar rhs_setup_data[2] = {R, l}; 184 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, 185 sizeof rhs_setup_data, &rhs_setup_data); 186 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 187 CeedQFunctionContextDestroy(&ctx_rhs_setup); 188 189 // Setup RHS and target 190 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 191 192 // Cleanup 193 CeedQFunctionDestroy(&qf_setup_rhs); 194 CeedOperatorDestroy(&op_setup_rhs); 195 } 196 197 // Cleanup 198 CeedQFunctionDestroy(&qf_setup_geo); 199 CeedOperatorDestroy(&op_setup_geo); 200 CeedVectorDestroy(&x_coord); 201 202 // Save libCEED data required for level 203 data->basis_x = basis_x; data->basis_u = basis_u; 204 data->elem_restr_x = elem_restr_x; 205 data->elem_restr_u = elem_restr_u; 206 data->elem_restr_u_i = elem_restr_u_i; 207 data->elem_restr_qd_i = elem_restr_qd_i; 208 data->qf_apply = qf_apply; 209 data->op_apply = op_apply; 210 data->q_data = q_data; 211 data->x_ceed = x_ceed; 212 data->y_ceed = y_ceed; 213 214 PetscFunctionReturn(0); 215 }; 216 217 // ----------------------------------------------------------------------------- 218 // Setup libCEED level transfer operator objects 219 // ----------------------------------------------------------------------------- 220 PetscErrorCode CeedLevelTransferSetup(DM dm, Ceed ceed, CeedInt level, 221 CeedInt num_comp_u, CeedData *data, 222 <<<<<<< HEAD 223 CeedInt *level_degrees, 224 CeedQFunction qf_restrict, CeedQFunction qf_prolong) { 225 PetscFunctionBeginUser; 226 // Return early if num_levels=1 227 if (num_levels == 1) 228 PetscFunctionReturn(0); 229 ======= 230 Vec fine_mult) { 231 int ierr; 232 >>>>>>> 2e9bde68 (Used "CeedOperatorMultigridLevelCreate" to create multigrid operators) 233 234 // Restriction - Fine to corse 235 CeedOperator op_restrict; 236 // Interpolation - Corse to fine 237 CeedOperator op_prolong; 238 // Coarse grid operator 239 CeedOperator op_apply; 240 // Basis 241 CeedBasis basis_u; 242 ierr = CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, &basis_u); 243 CHKERRQ(ierr); 244 245 // --------------------------------------------------------------------------- 246 // Coarse Grid, Prolongation, and Restriction Operators 247 // --------------------------------------------------------------------------- 248 // Create the Operators that compute the prolongation and 249 // restriction between the p-multigrid levels and the coarse grid eval. 250 // --------------------------------------------------------------------------- 251 // Place in libCEED array 252 const PetscScalar *m; 253 PetscMemType m_mem_type; 254 ierr = VecGetArrayReadAndMemType(fine_mult, &m, &m_mem_type); 255 CHKERRQ(ierr); 256 CeedVectorSetArray(data[level]->x_ceed, MemTypeP2C(m_mem_type), 257 CEED_USE_POINTER, (CeedScalar *)m); 258 259 CeedOperatorMultigridLevelCreate(data[level]->op_apply, data[level]->x_ceed, 260 data[level-1]->elem_restr_u, basis_u, 261 &op_apply, &op_prolong, &op_restrict); 262 263 // Restore PETSc vector 264 CeedVectorTakeArray(data[level]->x_ceed, MemTypeP2C(m_mem_type), 265 (CeedScalar **)&m); 266 ierr = VecRestoreArrayReadAndMemType(fine_mult, &m); CHKERRQ(ierr); 267 ierr = VecZeroEntries(fine_mult); CHKERRQ(ierr); 268 // -- Save libCEED data 269 data[level-1]->op_apply = op_apply; 270 data[level]->op_prolong = op_prolong; 271 data[level]->op_restrict = op_restrict; 272 273 CeedBasisDestroy(&basis_u); 274 PetscFunctionReturn(0); 275 }; 276 277 // ----------------------------------------------------------------------------- 278