1 #include "../include/swarmutils.h" 2 #include "../qfunctions/swarm/swarmmass.h" 3 4 // ------------------------------------------------------------------------------------------------ 5 // Context utilities 6 // ------------------------------------------------------------------------------------------------ 7 PetscErrorCode DMSwarmCeedContextCreate(DM dm_swarm, const char *ceed_resource, DMSwarmCeedContext *ctx) { 8 DM dm_mesh, dm_coord; 9 CeedElemRestriction elem_restr_u_mesh, elem_restr_x_mesh, elem_restr_x_points, elem_restr_u_points, elem_restr_q_data_points; 10 CeedBasis basis_u, basis_x; 11 CeedVector x_ref_points, q_data_points; 12 CeedInt num_comp; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscNew(ctx)); 16 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 17 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord)); 18 19 CeedInit(ceed_resource, &(*ctx)->ceed); 20 // Background mesh objects 21 { 22 BPData bp_data = {.q_mode = CEED_GAUSS}; 23 24 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &basis_u)); 25 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_coord, NULL, 0, 0, 0, bp_data, &basis_x)); 26 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_mesh, 0, NULL, 0, &elem_restr_u_mesh)); 27 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_coord, 0, NULL, 0, &elem_restr_x_mesh)); 28 29 // -- Mesh vectors 30 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->u_mesh, NULL); 31 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->v_mesh, NULL); 32 } 33 // Swarm objects 34 { 35 PetscInt dim; 36 const PetscInt *cell_points; 37 IS is_points; 38 Vec X_ref; 39 CeedInt num_elem; 40 41 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref)); 42 PetscCall(DMGetDimension(dm_mesh, &dim)); 43 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem); 44 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp); 45 46 PetscCall(ISGetIndices(is_points, &cell_points)); 47 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2; 48 CeedInt offsets[num_elem + 1 + num_points]; 49 50 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1; 51 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1]; 52 PetscCall(ISRestoreIndices(is_points, &cell_points)); 53 54 // -- Points restrictions 55 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 56 &elem_restr_u_points); 57 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 58 &elem_restr_x_points); 59 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 60 &elem_restr_q_data_points); 61 62 // -- Points vectors 63 CeedElemRestrictionCreateVector(elem_restr_u_points, &(*ctx)->u_points, NULL); 64 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL); 65 66 // -- Ref coordinates 67 { 68 PetscMemType X_mem_type; 69 const PetscScalar *x; 70 71 CeedVectorCreate((*ctx)->ceed, num_points * dim, &x_ref_points); 72 73 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type)); 74 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x); 75 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x)); 76 } 77 78 // Create Q data 79 { 80 CeedQFunction qf_setup; 81 CeedOperator op_setup; 82 CeedVector x_coord; 83 84 { 85 Vec X_loc; 86 CeedInt len; 87 const PetscScalar *x; 88 89 PetscCall(DMGetCoordinatesLocal(dm_mesh, &X_loc)); 90 PetscCall(VecGetLocalSize(X_loc, &len)); 91 CeedVectorCreate((*ctx)->ceed, len, &x_coord); 92 93 PetscCall(VecGetArrayRead(X_loc, &x)); 94 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)x); 95 PetscCall(VecRestoreArrayRead(X_loc, &x)); 96 } 97 98 // Setup geometric scaling 99 CeedQFunctionCreateInterior((*ctx)->ceed, 1, SetupMass, SetupMass_loc, &qf_setup); 100 CeedQFunctionAddInput(qf_setup, "x", dim * dim, CEED_EVAL_GRAD); 101 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 102 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 103 104 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 105 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 106 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 107 CeedOperatorSetField(op_setup, "rho", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 108 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points); 109 110 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE); 111 112 // Cleanup 113 CeedVectorDestroy(&x_coord); 114 CeedQFunctionDestroy(&qf_setup); 115 CeedOperatorDestroy(&op_setup); 116 } 117 118 // -- Cleanup 119 PetscCall(ISDestroy(&is_points)); 120 PetscCall(VecDestroy(&X_ref)); 121 } 122 123 PetscCall(DMSetApplicationContext(dm_mesh, (void *)(*ctx))); 124 125 // Create operators 126 // Mesh to points interpolation operator 127 { 128 CeedQFunction qf_mesh_to_points; 129 130 // -- Create operator 131 CeedQFunctionCreateIdentity((*ctx)->ceed, num_comp, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_mesh_to_points); 132 133 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mesh_to_points, NULL, NULL, &(*ctx)->op_mesh_to_points); 134 CeedOperatorSetField((*ctx)->op_mesh_to_points, "input", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 135 CeedOperatorSetField((*ctx)->op_mesh_to_points, "output", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 136 CeedOperatorAtPointsSetPoints((*ctx)->op_mesh_to_points, elem_restr_x_points, x_ref_points); 137 138 // -- Cleanup 139 CeedQFunctionDestroy(&qf_mesh_to_points); 140 } 141 142 // RHS operator 143 { 144 CeedQFunction qf_pts_to_mesh; 145 CeedQFunctionContext qf_ctx; 146 147 // -- Mass QFunction 148 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_pts_to_mesh); 149 CeedQFunctionAddInput(qf_pts_to_mesh, "q data", 1, CEED_EVAL_NONE); 150 CeedQFunctionAddInput(qf_pts_to_mesh, "u", num_comp, CEED_EVAL_NONE); 151 CeedQFunctionAddOutput(qf_pts_to_mesh, "v", num_comp, CEED_EVAL_INTERP); 152 153 // -- QFunction context 154 CeedQFunctionContextCreate((*ctx)->ceed, &qf_ctx); 155 CeedQFunctionContextSetData(qf_ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp); 156 CeedQFunctionSetContext(qf_pts_to_mesh, qf_ctx); 157 158 // -- Mass Operator 159 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_pts_to_mesh, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_points_to_mesh); 160 CeedOperatorSetField((*ctx)->op_points_to_mesh, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 161 CeedOperatorSetField((*ctx)->op_points_to_mesh, "u", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 162 CeedOperatorSetField((*ctx)->op_points_to_mesh, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 163 CeedOperatorAtPointsSetPoints((*ctx)->op_points_to_mesh, elem_restr_x_points, x_ref_points); 164 165 // -- Cleanup 166 CeedQFunctionContextDestroy(&qf_ctx); 167 CeedQFunctionDestroy(&qf_pts_to_mesh); 168 } 169 170 // Mass operator 171 { 172 CeedQFunction qf_mass; 173 CeedQFunctionContext ctx_mass; 174 175 // -- Mass QFunction 176 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_mass); 177 CeedQFunctionAddInput(qf_mass, "q data", 1, CEED_EVAL_NONE); 178 CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP); 179 CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP); 180 181 // -- QFunction context 182 CeedQFunctionContextCreate((*ctx)->ceed, &ctx_mass); 183 CeedQFunctionContextSetData(ctx_mass, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp); 184 CeedQFunctionSetContext(qf_mass, ctx_mass); 185 186 // -- Mass Operator 187 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_mass); 188 CeedOperatorSetField((*ctx)->op_mass, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 189 CeedOperatorSetField((*ctx)->op_mass, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 190 CeedOperatorSetField((*ctx)->op_mass, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 191 CeedOperatorAtPointsSetPoints((*ctx)->op_mass, elem_restr_x_points, x_ref_points); 192 193 // -- Cleanup 194 CeedQFunctionContextDestroy(&ctx_mass); 195 CeedQFunctionDestroy(&qf_mass); 196 } 197 198 // Cleanup 199 CeedElemRestrictionDestroy(&elem_restr_u_mesh); 200 CeedElemRestrictionDestroy(&elem_restr_x_mesh); 201 CeedElemRestrictionDestroy(&elem_restr_u_points); 202 CeedElemRestrictionDestroy(&elem_restr_x_points); 203 CeedElemRestrictionDestroy(&elem_restr_q_data_points); 204 CeedBasisDestroy(&basis_u); 205 CeedBasisDestroy(&basis_x); 206 CeedVectorDestroy(&x_ref_points); 207 CeedVectorDestroy(&q_data_points); 208 209 PetscFunctionReturn(PETSC_SUCCESS); 210 } 211 212 PetscErrorCode DMSwarmCeedContextDestroy(DMSwarmCeedContext *ctx) { 213 PetscFunctionBeginUser; 214 CeedDestroy(&(*ctx)->ceed); 215 CeedVectorDestroy(&(*ctx)->u_mesh); 216 CeedVectorDestroy(&(*ctx)->v_mesh); 217 CeedVectorDestroy(&(*ctx)->u_points); 218 CeedOperatorDestroy(&(*ctx)->op_mesh_to_points); 219 CeedOperatorDestroy(&(*ctx)->op_points_to_mesh); 220 CeedOperatorDestroy(&(*ctx)->op_mass); 221 PetscCall(PetscFree(*ctx)); 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 // ------------------------------------------------------------------------------------------------ 226 // PETSc-libCEED memory space utilities 227 // ------------------------------------------------------------------------------------------------ 228 PetscErrorCode DMSwarmPICFieldP2C(DM dm_swarm, const char *field, CeedVector x_ceed) { 229 PetscScalar *x; 230 231 PetscFunctionBeginUser; 232 PetscCall(DMSwarmGetField(dm_swarm, field, NULL, NULL, (void **)&x)); 233 CeedVectorSetArray(x_ceed, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)x); 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 PetscErrorCode DMSwarmPICFieldC2P(DM dm_swarm, const char *field, CeedVector x_ceed) { 238 PetscScalar *x; 239 240 PetscFunctionBeginUser; 241 CeedVectorTakeArray(x_ceed, CEED_MEM_HOST, (CeedScalar **)&x); 242 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&x)); 243 PetscFunctionReturn(PETSC_SUCCESS); 244 } 245 246 // ------------------------------------------------------------------------------------------------ 247 // Swarm point location utility 248 // ------------------------------------------------------------------------------------------------ 249 PetscErrorCode DMSwarmInitalizePointLocations(DM dm_swarm, PointSwarmType point_swarm_type, PetscInt num_points, PetscInt num_points_per_cell) { 250 PetscFunctionBeginUser; 251 252 switch (point_swarm_type) { 253 case SWARM_GAUSS: 254 case SWARM_UNIFORM: { 255 // -- Set gauss or uniform point locations in each cell 256 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)), dim = 3; 257 PetscScalar point_coords[num_points_per_cell * 3]; 258 CeedScalar points_1d[num_points_per_cell_1d], weights_1d[num_points_per_cell_1d]; 259 260 if (point_swarm_type == SWARM_GAUSS) { 261 PetscCall(CeedGaussQuadrature(num_points_per_cell_1d, points_1d, weights_1d)); 262 } else { 263 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) points_1d[i] = 2.0 * (PetscReal)(i + 1) / (PetscReal)(num_points_per_cell_1d + 1) - 1; 264 } 265 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) { 266 for (PetscInt j = 0; j < num_points_per_cell_1d; j++) { 267 for (PetscInt k = 0; k < num_points_per_cell_1d; k++) { 268 PetscInt p = (i * num_points_per_cell_1d + j) * num_points_per_cell_1d + k; 269 270 point_coords[p * dim + 0] = points_1d[i]; 271 point_coords[p * dim + 1] = points_1d[j]; 272 point_coords[p * dim + 2] = points_1d[k]; 273 } 274 } 275 } 276 PetscCall(DMSwarmSetPointCoordinatesCellwise(dm_swarm, num_points_per_cell_1d * num_points_per_cell_1d * num_points_per_cell_1d, point_coords)); 277 } break; 278 case SWARM_CELL_RANDOM: { 279 // -- Set points randomly in each cell 280 PetscInt dim = 3, num_cells_total = 1, num_cells[] = {1, 1, 1}; 281 PetscScalar *point_coords; 282 PetscRandom rng; 283 284 PetscOptionsBegin(PetscObjectComm((PetscObject)dm_swarm), NULL, "libCEED example using PETSc with DMSwarm", NULL); 285 286 PetscCall(PetscOptionsInt("-dm_plex_dim", "Background mesh dimension", NULL, dim, &dim, NULL)); 287 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of cells", NULL, num_cells, &dim, NULL)); 288 289 PetscOptionsEnd(); 290 291 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm_swarm), &rng)); 292 293 num_cells_total = num_cells[0] * num_cells[1] * num_cells[2]; 294 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 295 for (PetscInt c = 0; c < num_cells_total; c++) { 296 PetscInt cell_index[3] = {c % num_cells[0], (c / num_cells[0]) % num_cells[1], (c / num_cells[0] / num_cells[1]) % num_cells[2]}; 297 298 for (PetscInt p = 0; p < num_points_per_cell; p++) { 299 PetscInt point_index = c * num_points_per_cell + p; 300 PetscScalar random_value; 301 302 for (PetscInt i = 0; i < dim; i++) { 303 PetscCall(PetscRandomGetValue(rng, &random_value)); 304 point_coords[point_index * dim + i] = -1.0 + cell_index[i] * 2.0 / (num_cells[i] + 1.0) + random_value; 305 } 306 } 307 } 308 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 309 PetscCall(PetscRandomDestroy(&rng)); 310 } break; 311 case SWARM_SINUSOIDAL: { 312 // -- Set points distributed per sinusoidal functions 313 PetscInt dim = 3; 314 PetscScalar *point_coords; 315 DM dm_mesh; 316 317 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 318 PetscCall(DMGetDimension(dm_mesh, &dim)); 319 320 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 321 for (PetscInt p = 0; p < num_points; p++) { 322 point_coords[p * dim + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 323 if (dim > 1) point_coords[p * dim + 1] = -PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 324 if (dim > 2) point_coords[p * dim + 2] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 325 } 326 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 327 } break; 328 } 329 PetscCall(DMSwarmMigrate(dm_swarm, PETSC_TRUE)); 330 PetscFunctionReturn(PETSC_SUCCESS); 331 } 332 333 /*@C 334 DMSwarmCreateReferenceCoordinates - Compute the cell reference coordinates for local DMSwarm points. 335 336 Collective 337 338 Input Parameter: 339 . dm_swarm - the `DMSwarm` 340 341 Output Parameters: 342 + is_points - The IS object for indexing into points per cell 343 - X_points_ref - Vec holding the cell reference coordinates for local DMSwarm points 344 345 The index set contains ranges of indices for each local cell. This range contains the indices of every point in the cell. 346 347 ``` 348 total_num_cells 349 cell_0_start_index 350 cell_1_start_index 351 cell_2_start_index 352 ... 353 cell_n_start_index 354 cell_n_stop_index 355 cell_0_point_0 356 cell_0_point_0 357 ... 358 cell_n_point_m 359 ``` 360 361 Level: beginner 362 363 .seealso: `DMSwarm` 364 @*/ 365 PetscErrorCode DMSwarmCreateReferenceCoordinates(DM dm_swarm, IS *is_points, Vec *X_points_ref) { 366 PetscInt cell_start, cell_end, num_cells_local, dim, num_points_local, *cell_points, points_offset; 367 PetscScalar *coords_points_ref; 368 const PetscScalar *coords_points_true; 369 DM dm_mesh; 370 371 PetscFunctionBeginUser; 372 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 373 374 // Create vector to hold reference coordinates 375 { 376 Vec X_points_true; 377 378 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 379 PetscCall(VecDuplicate(X_points_true, X_points_ref)); 380 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 381 } 382 383 // Allocate index set array 384 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 385 num_cells_local = cell_end - cell_start; 386 points_offset = num_cells_local + 2; 387 PetscCall(VecGetLocalSize(*X_points_ref, &num_points_local)); 388 PetscCall(DMGetDimension(dm_mesh, &dim)); 389 num_points_local /= dim; 390 PetscCall(PetscMalloc1(num_points_local + num_cells_local + 2, &cell_points)); 391 cell_points[0] = num_cells_local; 392 393 // Get reference coordinates for each swarm point wrt the elements in the background mesh 394 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 395 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 396 PetscCall(VecGetArray(*X_points_ref, &coords_points_ref)); 397 for (PetscInt cell = cell_start, num_points_processed = 0; cell < cell_end; cell++) { 398 PetscInt *points_in_cell, num_points_in_cell, local_cell = cell - cell_start; 399 PetscReal v[3], J[9], invJ[9], detJ, v0_ref[3] = {-1.0, -1.0, -1.0}; 400 401 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell)); 402 // -- Reference coordinates for swarm points in background mesh element 403 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ)); 404 cell_points[local_cell + 1] = num_points_processed + points_offset; 405 for (PetscInt p = 0; p < num_points_in_cell; p++) { 406 const CeedInt point = points_in_cell[p]; 407 408 cell_points[points_offset + (num_points_processed++)] = point; 409 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_true[point * dim], &coords_points_ref[point * dim]); 410 } 411 412 // -- Cleanup 413 PetscCall(PetscFree(points_in_cell)); 414 } 415 cell_points[points_offset - 1] = num_points_local + points_offset; 416 417 // Cleanup 418 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 419 PetscCall(VecRestoreArray(*X_points_ref, &coords_points_ref)); 420 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 421 422 // Create index set 423 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_points_local + points_offset, cell_points, PETSC_OWN_POINTER, is_points)); 424 PetscFunctionReturn(PETSC_SUCCESS); 425 } 426 427 // ------------------------------------------------------------------------------------------------ 428 // RHS for Swarm to Mesh projection 429 // ------------------------------------------------------------------------------------------------ 430 PetscErrorCode DMSwarmCreateProjectionRHS(DM dm_swarm, const char *field, Vec B_mesh) { 431 PetscMemType B_mem_type; 432 DM dm_mesh; 433 Vec B_mesh_loc; 434 DMSwarmCeedContext swarm_ceed_context; 435 436 PetscFunctionBeginUser; 437 // Get mesh DM 438 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 439 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 440 441 // Get mesh values 442 PetscCall(DMGetLocalVector(dm_mesh, &B_mesh_loc)); 443 PetscCall(VecZeroEntries(B_mesh_loc)); 444 PetscCall(VecP2C(B_mesh_loc, &B_mem_type, swarm_ceed_context->v_mesh)); 445 446 // Get swarm access 447 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 448 PetscCall(DMSwarmPICFieldP2C(dm_swarm, field, swarm_ceed_context->u_points)); 449 450 // Interpolate field from swarm points to mesh 451 CeedOperatorApply(swarm_ceed_context->op_points_to_mesh, swarm_ceed_context->u_points, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE); 452 453 // Restore PETSc Vecs and Local to Global 454 PetscCall(VecC2P(swarm_ceed_context->v_mesh, B_mem_type, B_mesh_loc)); 455 PetscCall(VecZeroEntries(B_mesh)); 456 PetscCall(DMLocalToGlobal(dm_mesh, B_mesh_loc, ADD_VALUES, B_mesh)); 457 458 // Cleanup 459 PetscCall(DMSwarmPICFieldC2P(dm_swarm, field, swarm_ceed_context->u_points)); 460 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 461 PetscCall(DMRestoreLocalVector(dm_mesh, &B_mesh_loc)); 462 PetscFunctionReturn(PETSC_SUCCESS); 463 } 464 465 // ------------------------------------------------------------------------------------------------ 466 // Swarm "mass matrix" 467 // ------------------------------------------------------------------------------------------------ 468 PetscErrorCode MatMult_SwarmMass(Mat A, Vec U_mesh, Vec V_mesh) { 469 PetscMemType U_mem_type, V_mem_type; 470 DM dm_mesh; 471 Vec U_mesh_loc, V_mesh_loc; 472 DMSwarmCeedContext swarm_ceed_context; 473 474 PetscFunctionBeginUser; 475 // Get mesh DM 476 PetscCall(MatGetDM(A, &dm_mesh)); 477 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 478 479 // Global to Local and get PETSc Vec access 480 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc)); 481 PetscCall(VecZeroEntries(U_mesh_loc)); 482 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc)); 483 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh)); 484 PetscCall(DMGetLocalVector(dm_mesh, &V_mesh_loc)); 485 PetscCall(VecZeroEntries(V_mesh_loc)); 486 PetscCall(VecP2C(V_mesh_loc, &V_mem_type, swarm_ceed_context->v_mesh)); 487 488 // Apply swarm mass operator 489 CeedOperatorApply(swarm_ceed_context->op_mass, swarm_ceed_context->u_mesh, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE); 490 491 // Restore PETSc Vecs and Local to Global 492 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh, U_mem_type, U_mesh_loc)); 493 PetscCall(VecC2P(swarm_ceed_context->v_mesh, V_mem_type, V_mesh_loc)); 494 PetscCall(VecZeroEntries(V_mesh)); 495 PetscCall(DMLocalToGlobal(dm_mesh, V_mesh_loc, ADD_VALUES, V_mesh)); 496 497 // Cleanup 498 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc)); 499 PetscCall(DMRestoreLocalVector(dm_mesh, &V_mesh_loc)); 500 PetscFunctionReturn(PETSC_SUCCESS); 501 } 502 503 // ------------------------------------------------------------------------------------------------ 504 // Swarm to mesh projection 505 // ------------------------------------------------------------------------------------------------ 506 PetscErrorCode DMSwarmProjectFromSwarmToCells(DM dm_swarm, const char *field, Vec U_mesh) { 507 PetscBool test_mode; 508 Vec B_mesh; 509 Mat M; 510 KSP ksp; 511 DM dm_mesh; 512 DMSwarmCeedContext swarm_ceed_context; 513 MPI_Comm comm; 514 515 PetscFunctionBeginUser; 516 517 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swarm-to-Mesh Projection Options", NULL); 518 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 519 PetscOptionsEnd(); 520 521 comm = PetscObjectComm((PetscObject)dm_swarm); 522 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 523 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 524 PetscCall(VecDuplicate(U_mesh, &B_mesh)); 525 526 // Setup "mass matrix" 527 { 528 PetscInt l_size, g_size; 529 530 PetscCall(VecGetLocalSize(U_mesh, &l_size)); 531 PetscCall(VecGetSize(U_mesh, &g_size)); 532 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, swarm_ceed_context, &M)); 533 PetscCall(MatSetDM(M, dm_mesh)); 534 PetscCall(MatShellSetOperation(M, MATOP_MULT, (void (*)(void))MatMult_SwarmMass)); 535 } 536 537 // Setup KSP 538 { 539 PC pc; 540 541 PetscCall(KSPCreate(comm, &ksp)); 542 PetscCall(KSPGetPC(ksp, &pc)); 543 PetscCall(PCSetType(pc, PCJACOBI)); 544 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 545 PetscCall(KSPSetType(ksp, KSPCG)); 546 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 547 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 548 PetscCall(KSPSetOperators(ksp, M, M)); 549 PetscCall(KSPSetFromOptions(ksp)); 550 PetscCall(PetscObjectSetName((PetscObject)ksp, "Swarm-to-Mesh Projection")); 551 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_projection_view")); 552 } 553 554 // Setup RHS 555 PetscCall(DMSwarmCreateProjectionRHS(dm_swarm, field, B_mesh)); 556 557 // Solve 558 PetscCall(VecZeroEntries(U_mesh)); 559 PetscCall(KSPSolve(ksp, B_mesh, U_mesh)); 560 561 // KSP summary 562 { 563 KSPType ksp_type; 564 KSPConvergedReason reason; 565 PetscReal rnorm; 566 PetscInt its; 567 PetscCall(KSPGetType(ksp, &ksp_type)); 568 PetscCall(KSPGetConvergedReason(ksp, &reason)); 569 PetscCall(KSPGetIterationNumber(ksp, &its)); 570 PetscCall(KSPGetResidualNorm(ksp, &rnorm)); 571 572 if (!test_mode || reason < 0 || rnorm > 1e-8) { 573 PetscCall(PetscPrintf(comm, 574 "Swarm-to-Mesh Projection KSP Solve:\n" 575 " KSP type: %s\n" 576 " KSP convergence: %s\n" 577 " Total KSP iterations: %" PetscInt_FMT "\n" 578 " Final rnorm: %e\n", 579 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm)); 580 } 581 } 582 583 // Optional viewing 584 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_view")); 585 586 // Cleanup 587 PetscCall(VecDestroy(&B_mesh)); 588 PetscCall(MatDestroy(&M)); 589 PetscCall(KSPDestroy(&ksp)); 590 591 PetscFunctionReturn(PETSC_SUCCESS); 592 } 593