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