1 // Copyright (c) 2017-2024, 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 "../include/matops.h" 10 #include "../qfunctions/swarm/swarmmass.h" 11 12 // ------------------------------------------------------------------------------------------------ 13 // Context utilities 14 // ------------------------------------------------------------------------------------------------ 15 PetscErrorCode DMSwarmCeedContextCreate(DM dm_swarm, const char *ceed_resource, DMSwarmCeedContext *ctx) { 16 DM dm_mesh, dm_coord; 17 CeedElemRestriction elem_restr_u_mesh, elem_restr_x_mesh, elem_restr_x_points, elem_restr_u_points, elem_restr_q_data_points; 18 CeedBasis basis_u, basis_x; 19 CeedVector x_ref_points, q_data_points; 20 CeedInt num_comp; 21 22 PetscFunctionBeginUser; 23 PetscCall(PetscNew(ctx)); 24 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 25 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord)); 26 27 CeedInit(ceed_resource, &(*ctx)->ceed); 28 // Background mesh objects 29 { 30 BPData bp_data = {.q_mode = CEED_GAUSS}; 31 32 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &basis_u)); 33 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_coord, NULL, 0, 0, 0, bp_data, &basis_x)); 34 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_mesh, 0, NULL, 0, &elem_restr_u_mesh)); 35 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_coord, 0, NULL, 0, &elem_restr_x_mesh)); 36 37 // -- Mesh vectors 38 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->u_mesh, NULL); 39 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->v_mesh, NULL); 40 } 41 // Swarm objects 42 { 43 PetscInt dim; 44 const PetscInt *cell_points; 45 IS is_points; 46 Vec X_ref; 47 CeedInt num_elem; 48 49 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref)); 50 PetscCall(DMGetDimension(dm_mesh, &dim)); 51 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem); 52 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp); 53 54 PetscCall(ISGetIndices(is_points, &cell_points)); 55 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2; 56 CeedInt offsets[num_elem + 1 + num_points]; 57 58 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1; 59 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1]; 60 PetscCall(ISRestoreIndices(is_points, &cell_points)); 61 62 // -- Points restrictions 63 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 64 &elem_restr_u_points); 65 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 66 &elem_restr_x_points); 67 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 68 &elem_restr_q_data_points); 69 70 // -- Points vectors 71 CeedElemRestrictionCreateVector(elem_restr_u_points, &(*ctx)->u_points, NULL); 72 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL); 73 74 // -- Ref coordinates 75 { 76 PetscMemType X_mem_type; 77 const PetscScalar *x; 78 79 CeedVectorCreate((*ctx)->ceed, num_points * dim, &x_ref_points); 80 81 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type)); 82 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x); 83 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x)); 84 } 85 86 // Create Q data 87 { 88 CeedQFunction qf_setup; 89 CeedOperator op_setup; 90 CeedVector x_coord; 91 92 { 93 Vec X_loc; 94 PetscInt len; 95 const PetscScalar *x; 96 97 PetscCall(DMGetCoordinatesLocal(dm_mesh, &X_loc)); 98 PetscCall(VecGetLocalSize(X_loc, &len)); 99 CeedVectorCreate((*ctx)->ceed, len, &x_coord); 100 101 PetscCall(VecGetArrayRead(X_loc, &x)); 102 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)x); 103 PetscCall(VecRestoreArrayRead(X_loc, &x)); 104 } 105 106 // Setup geometric scaling 107 CeedQFunctionCreateInterior((*ctx)->ceed, 1, SetupMass, SetupMass_loc, &qf_setup); 108 CeedQFunctionAddInput(qf_setup, "x", dim * dim, CEED_EVAL_GRAD); 109 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 110 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE); 111 112 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 113 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 114 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 115 CeedOperatorSetField(op_setup, "rho", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 116 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points); 117 118 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE); 119 120 // Cleanup 121 CeedVectorDestroy(&x_coord); 122 CeedQFunctionDestroy(&qf_setup); 123 CeedOperatorDestroy(&op_setup); 124 } 125 126 // -- Cleanup 127 PetscCall(ISDestroy(&is_points)); 128 PetscCall(VecDestroy(&X_ref)); 129 } 130 131 PetscCall(DMSetApplicationContext(dm_mesh, (void *)(*ctx))); 132 133 // Create operators 134 // Mesh to points interpolation operator 135 { 136 CeedQFunction qf_mesh_to_points; 137 138 // -- Create operator 139 CeedQFunctionCreateIdentity((*ctx)->ceed, num_comp, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_mesh_to_points); 140 141 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mesh_to_points, NULL, NULL, &(*ctx)->op_mesh_to_points); 142 CeedOperatorSetField((*ctx)->op_mesh_to_points, "input", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 143 CeedOperatorSetField((*ctx)->op_mesh_to_points, "output", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 144 CeedOperatorAtPointsSetPoints((*ctx)->op_mesh_to_points, elem_restr_x_points, x_ref_points); 145 146 // -- Cleanup 147 CeedQFunctionDestroy(&qf_mesh_to_points); 148 } 149 150 // RHS operator 151 { 152 CeedQFunction qf_pts_to_mesh; 153 CeedQFunctionContext qf_ctx; 154 155 // -- Mass QFunction 156 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_pts_to_mesh); 157 CeedQFunctionAddInput(qf_pts_to_mesh, "q data", 1, CEED_EVAL_NONE); 158 CeedQFunctionAddInput(qf_pts_to_mesh, "u", num_comp, CEED_EVAL_NONE); 159 CeedQFunctionAddOutput(qf_pts_to_mesh, "v", num_comp, CEED_EVAL_INTERP); 160 161 // -- QFunction context 162 CeedQFunctionContextCreate((*ctx)->ceed, &qf_ctx); 163 CeedQFunctionContextSetData(qf_ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp); 164 CeedQFunctionSetContext(qf_pts_to_mesh, qf_ctx); 165 166 // -- Mass Operator 167 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_pts_to_mesh, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_points_to_mesh); 168 CeedOperatorSetField((*ctx)->op_points_to_mesh, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 169 CeedOperatorSetField((*ctx)->op_points_to_mesh, "u", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 170 CeedOperatorSetField((*ctx)->op_points_to_mesh, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 171 CeedOperatorAtPointsSetPoints((*ctx)->op_points_to_mesh, elem_restr_x_points, x_ref_points); 172 173 // -- Cleanup 174 CeedQFunctionContextDestroy(&qf_ctx); 175 CeedQFunctionDestroy(&qf_pts_to_mesh); 176 } 177 178 // Mass operator 179 { 180 CeedQFunction qf_mass; 181 CeedQFunctionContext ctx_mass; 182 183 // -- Mass QFunction 184 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_mass); 185 CeedQFunctionAddInput(qf_mass, "q data", 1, CEED_EVAL_NONE); 186 CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP); 187 CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP); 188 189 // -- QFunction context 190 CeedQFunctionContextCreate((*ctx)->ceed, &ctx_mass); 191 CeedQFunctionContextSetData(ctx_mass, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp); 192 CeedQFunctionSetContext(qf_mass, ctx_mass); 193 194 // -- Mass Operator 195 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_mass); 196 CeedOperatorSetField((*ctx)->op_mass, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 197 CeedOperatorSetField((*ctx)->op_mass, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 198 CeedOperatorSetField((*ctx)->op_mass, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 199 CeedOperatorAtPointsSetPoints((*ctx)->op_mass, elem_restr_x_points, x_ref_points); 200 201 // -- Cleanup 202 CeedQFunctionContextDestroy(&ctx_mass); 203 CeedQFunctionDestroy(&qf_mass); 204 } 205 206 // Cleanup 207 CeedElemRestrictionDestroy(&elem_restr_u_mesh); 208 CeedElemRestrictionDestroy(&elem_restr_x_mesh); 209 CeedElemRestrictionDestroy(&elem_restr_u_points); 210 CeedElemRestrictionDestroy(&elem_restr_x_points); 211 CeedElemRestrictionDestroy(&elem_restr_q_data_points); 212 CeedBasisDestroy(&basis_u); 213 CeedBasisDestroy(&basis_x); 214 CeedVectorDestroy(&x_ref_points); 215 CeedVectorDestroy(&q_data_points); 216 PetscFunctionReturn(PETSC_SUCCESS); 217 } 218 219 PetscErrorCode DMSwarmCeedContextDestroy(DMSwarmCeedContext *ctx) { 220 PetscFunctionBeginUser; 221 CeedDestroy(&(*ctx)->ceed); 222 CeedVectorDestroy(&(*ctx)->u_mesh); 223 CeedVectorDestroy(&(*ctx)->v_mesh); 224 CeedVectorDestroy(&(*ctx)->u_points); 225 CeedOperatorDestroy(&(*ctx)->op_mesh_to_points); 226 CeedOperatorDestroy(&(*ctx)->op_points_to_mesh); 227 CeedOperatorDestroy(&(*ctx)->op_mass); 228 PetscCall(PetscFree(*ctx)); 229 PetscFunctionReturn(PETSC_SUCCESS); 230 } 231 232 // ------------------------------------------------------------------------------------------------ 233 // PETSc-libCEED memory space utilities 234 // ------------------------------------------------------------------------------------------------ 235 PetscErrorCode DMSwarmPICFieldP2C(DM dm_swarm, const char *field, CeedVector x_ceed) { 236 PetscScalar *x; 237 238 PetscFunctionBeginUser; 239 PetscCall(DMSwarmGetField(dm_swarm, field, NULL, NULL, (void **)&x)); 240 CeedVectorSetArray(x_ceed, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)x); 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 PetscErrorCode DMSwarmPICFieldC2P(DM dm_swarm, const char *field, CeedVector x_ceed) { 245 PetscScalar *x; 246 247 PetscFunctionBeginUser; 248 CeedVectorTakeArray(x_ceed, CEED_MEM_HOST, (CeedScalar **)&x); 249 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&x)); 250 PetscFunctionReturn(PETSC_SUCCESS); 251 } 252 253 // ------------------------------------------------------------------------------------------------ 254 // Swarm point location utility 255 // ------------------------------------------------------------------------------------------------ 256 PetscErrorCode DMSwarmInitalizePointLocations(DM dm_swarm, PointSwarmType point_swarm_type, PetscInt num_points, PetscInt num_points_per_cell) { 257 PetscFunctionBeginUser; 258 switch (point_swarm_type) { 259 case SWARM_GAUSS: 260 case SWARM_UNIFORM: { 261 // -- Set gauss or uniform point locations in each cell 262 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)), dim = 3; 263 PetscScalar point_coords[num_points_per_cell * 3]; 264 CeedScalar points_1d[num_points_per_cell_1d], weights_1d[num_points_per_cell_1d]; 265 266 if (point_swarm_type == SWARM_GAUSS) { 267 PetscCall(CeedGaussQuadrature(num_points_per_cell_1d, points_1d, weights_1d)); 268 } else { 269 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) points_1d[i] = 2.0 * (PetscReal)(i + 0.5) / (PetscReal)num_points_per_cell_1d - 1; 270 } 271 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) { 272 for (PetscInt j = 0; j < num_points_per_cell_1d; j++) { 273 for (PetscInt k = 0; k < num_points_per_cell_1d; k++) { 274 PetscInt p = (i * num_points_per_cell_1d + j) * num_points_per_cell_1d + k; 275 276 point_coords[p * dim + 0] = points_1d[i]; 277 point_coords[p * dim + 1] = points_1d[j]; 278 point_coords[p * dim + 2] = points_1d[k]; 279 } 280 } 281 } 282 PetscCall(DMSwarmSetPointCoordinatesCellwise(dm_swarm, num_points_per_cell_1d * num_points_per_cell_1d * num_points_per_cell_1d, point_coords)); 283 } break; 284 case SWARM_CELL_RANDOM: { 285 DM dm_mesh; 286 287 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 288 // DMSwarmSetPointCoordinatesRandom expects local coordinates to be set up to ensure call is non-collective 289 PetscCall(DMGetCoordinatesLocalSetUp(dm_mesh)); 290 PetscCall(DMSwarmSetPointCoordinatesRandom(dm_swarm, num_points_per_cell)); 291 } break; 292 case SWARM_SINUSOIDAL: { 293 // -- Set points distributed per sinusoidal functions 294 PetscInt dim = 3; 295 PetscScalar *point_coords; 296 DM dm_mesh; 297 298 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 299 PetscCall(DMGetDimension(dm_mesh, &dim)); 300 301 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 302 for (PetscInt p = 0; p < num_points; p++) { 303 point_coords[p * dim + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 304 if (dim > 1) point_coords[p * dim + 1] = -PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 305 if (dim > 2) point_coords[p * dim + 2] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 306 } 307 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 308 } break; 309 } 310 PetscCall(DMSwarmMigrate(dm_swarm, PETSC_TRUE)); 311 PetscFunctionReturn(PETSC_SUCCESS); 312 } 313 314 /*@C 315 DMSwarmCreateReferenceCoordinates - Compute the cell reference coordinates for local DMSwarm points. 316 317 Collective 318 319 Input Parameter: 320 . dm_swarm - the `DMSwarm` 321 322 Output Parameters: 323 + is_points - The IS object for indexing into points per cell 324 - X_points_ref - Vec holding the cell reference coordinates for local DMSwarm points 325 326 The index set contains ranges of indices for each local cell. This range contains the indices of every point in the cell. 327 328 ``` 329 total_num_cells 330 cell_0_start_index 331 cell_1_start_index 332 cell_2_start_index 333 ... 334 cell_n_start_index 335 cell_n_stop_index 336 cell_0_point_0 337 cell_0_point_0 338 ... 339 cell_n_point_m 340 ``` 341 342 Level: beginner 343 344 .seealso: `DMSwarm` 345 @*/ 346 PetscErrorCode DMSwarmCreateReferenceCoordinates(DM dm_swarm, IS *is_points, Vec *X_points_ref) { 347 PetscInt cell_start, cell_end, num_cells_local, dim, num_points_local, *cell_points, points_offset; 348 PetscScalar *coords_points_ref; 349 const PetscScalar *coords_points_true; 350 DM dm_mesh; 351 352 PetscFunctionBeginUser; 353 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 354 355 // Create vector to hold reference coordinates 356 { 357 Vec X_points_true; 358 359 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 360 PetscCall(VecDuplicate(X_points_true, X_points_ref)); 361 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 362 } 363 364 // Allocate index set array 365 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 366 num_cells_local = cell_end - cell_start; 367 points_offset = num_cells_local + 2; 368 PetscCall(VecGetLocalSize(*X_points_ref, &num_points_local)); 369 PetscCall(DMGetDimension(dm_mesh, &dim)); 370 num_points_local /= dim; 371 PetscCall(PetscMalloc1(num_points_local + num_cells_local + 2, &cell_points)); 372 cell_points[0] = num_cells_local; 373 374 // Get reference coordinates for each swarm point wrt the elements in the background mesh 375 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 376 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 377 PetscCall(VecGetArray(*X_points_ref, &coords_points_ref)); 378 for (PetscInt cell = cell_start, num_points_processed = 0; cell < cell_end; cell++) { 379 PetscInt *points_in_cell, num_points_in_cell, local_cell = cell - cell_start; 380 PetscReal v[3], J[9], invJ[9], detJ, v0_ref[3] = {-1.0, -1.0, -1.0}; 381 382 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell)); 383 // -- Reference coordinates for swarm points in background mesh element 384 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ)); 385 cell_points[local_cell + 1] = num_points_processed + points_offset; 386 for (PetscInt p = 0; p < num_points_in_cell; p++) { 387 const CeedInt point = points_in_cell[p]; 388 389 cell_points[points_offset + (num_points_processed++)] = point; 390 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_true[point * dim], &coords_points_ref[point * dim]); 391 } 392 393 // -- Cleanup 394 PetscCall(DMSwarmSortRestorePointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell)); 395 } 396 cell_points[points_offset - 1] = num_points_local + points_offset; 397 398 // Cleanup 399 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 400 PetscCall(VecRestoreArray(*X_points_ref, &coords_points_ref)); 401 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 402 403 // Create index set 404 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_points_local + points_offset, cell_points, PETSC_OWN_POINTER, is_points)); 405 PetscFunctionReturn(PETSC_SUCCESS); 406 } 407 408 // ------------------------------------------------------------------------------------------------ 409 // RHS for Swarm to Mesh projection 410 // ------------------------------------------------------------------------------------------------ 411 PetscErrorCode DMSwarmCreateProjectionRHS(DM dm_swarm, const char *field, Vec U_points, Vec B_mesh) { 412 PetscMemType B_mem_type, U_mem_type; 413 DM dm_mesh; 414 Vec B_mesh_loc; 415 PetscBool has_u_points; 416 DMSwarmCeedContext swarm_ceed_context; 417 418 PetscFunctionBeginUser; 419 // Get mesh DM 420 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 421 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 422 423 // Get swarm access 424 has_u_points = !!U_points; 425 if (!has_u_points) { 426 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 427 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, field, &U_points)); 428 } 429 430 // Get mesh values 431 PetscCall(VecReadP2C(U_points, &U_mem_type, swarm_ceed_context->u_points)); 432 PetscCall(DMGetLocalVector(dm_mesh, &B_mesh_loc)); 433 PetscCall(VecZeroEntries(B_mesh_loc)); 434 PetscCall(VecP2C(B_mesh_loc, &B_mem_type, swarm_ceed_context->v_mesh)); 435 436 // Interpolate field from swarm points to mesh 437 CeedOperatorApply(swarm_ceed_context->op_points_to_mesh, swarm_ceed_context->u_points, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE); 438 439 // Restore PETSc Vecs and Local to Global 440 PetscCall(VecReadC2P(swarm_ceed_context->u_points, U_mem_type, U_points)); 441 PetscCall(VecC2P(swarm_ceed_context->v_mesh, B_mem_type, B_mesh_loc)); 442 PetscCall(VecZeroEntries(B_mesh)); 443 PetscCall(DMLocalToGlobal(dm_mesh, B_mesh_loc, ADD_VALUES, B_mesh)); 444 445 // Restore swarm access 446 if (!has_u_points) { 447 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, field, &U_points)); 448 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 449 } 450 451 // Cleanup 452 PetscCall(DMRestoreLocalVector(dm_mesh, &B_mesh_loc)); 453 PetscFunctionReturn(PETSC_SUCCESS); 454 } 455 456 // ------------------------------------------------------------------------------------------------ 457 // Swarm "mass matrix" 458 // ------------------------------------------------------------------------------------------------ 459 PetscErrorCode MatMult_SwarmMass(Mat A, Vec U_mesh, Vec V_mesh) { 460 PetscMemType U_mem_type, V_mem_type; 461 DM dm_mesh; 462 Vec U_mesh_loc, V_mesh_loc; 463 DMSwarmCeedContext swarm_ceed_context; 464 465 PetscFunctionBeginUser; 466 // Get mesh DM 467 PetscCall(MatGetDM(A, &dm_mesh)); 468 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 469 470 // Global to Local and get PETSc Vec access 471 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc)); 472 PetscCall(VecZeroEntries(U_mesh_loc)); 473 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc)); 474 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh)); 475 PetscCall(DMGetLocalVector(dm_mesh, &V_mesh_loc)); 476 PetscCall(VecZeroEntries(V_mesh_loc)); 477 PetscCall(VecP2C(V_mesh_loc, &V_mem_type, swarm_ceed_context->v_mesh)); 478 479 // Apply swarm mass operator 480 CeedOperatorApply(swarm_ceed_context->op_mass, swarm_ceed_context->u_mesh, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE); 481 482 // Restore PETSc Vecs and Local to Global 483 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh, U_mem_type, U_mesh_loc)); 484 PetscCall(VecC2P(swarm_ceed_context->v_mesh, V_mem_type, V_mesh_loc)); 485 PetscCall(VecZeroEntries(V_mesh)); 486 PetscCall(DMLocalToGlobal(dm_mesh, V_mesh_loc, ADD_VALUES, V_mesh)); 487 488 // Cleanup 489 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc)); 490 PetscCall(DMRestoreLocalVector(dm_mesh, &V_mesh_loc)); 491 PetscFunctionReturn(PETSC_SUCCESS); 492 } 493 494 // ------------------------------------------------------------------------------------------------ 495 // Swarm to mesh projection 496 // ------------------------------------------------------------------------------------------------ 497 PetscErrorCode DMSwarmProjectFromSwarmToCells(DM dm_swarm, const char *field, Vec U_points, Vec U_mesh) { 498 PetscBool test_mode; 499 Vec B_mesh; 500 Mat M; 501 KSP ksp; 502 DM dm_mesh; 503 DMSwarmCeedContext swarm_ceed_context; 504 MPI_Comm comm; 505 506 PetscFunctionBeginUser; 507 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swarm-to-Mesh Projection Options", NULL); 508 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 509 PetscOptionsEnd(); 510 511 comm = PetscObjectComm((PetscObject)dm_swarm); 512 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 513 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 514 PetscCall(VecDuplicate(U_mesh, &B_mesh)); 515 516 // Setup "mass matrix" 517 { 518 PetscInt l_size, g_size; 519 520 PetscCall(VecGetLocalSize(U_mesh, &l_size)); 521 PetscCall(VecGetSize(U_mesh, &g_size)); 522 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, swarm_ceed_context, &M)); 523 PetscCall(MatSetDM(M, dm_mesh)); 524 PetscCall(MatShellSetOperation(M, MATOP_MULT, (void (*)(void))MatMult_SwarmMass)); 525 } 526 527 // Setup KSP 528 { 529 PC pc; 530 531 PetscCall(KSPCreate(comm, &ksp)); 532 PetscCall(KSPGetPC(ksp, &pc)); 533 PetscCall(PCSetType(pc, PCJACOBI)); 534 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 535 PetscCall(KSPSetType(ksp, KSPCG)); 536 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 537 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 538 PetscCall(KSPSetOperators(ksp, M, M)); 539 PetscCall(KSPSetFromOptions(ksp)); 540 PetscCall(PetscObjectSetName((PetscObject)ksp, "Swarm-to-Mesh Projection")); 541 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_projection_view")); 542 } 543 544 // Setup RHS 545 PetscCall(DMSwarmCreateProjectionRHS(dm_swarm, field, U_points, B_mesh)); 546 547 // Solve 548 PetscCall(VecZeroEntries(U_mesh)); 549 PetscCall(KSPSolve(ksp, B_mesh, U_mesh)); 550 551 // KSP summary 552 { 553 KSPType ksp_type; 554 KSPConvergedReason reason; 555 PetscReal rnorm; 556 PetscInt its; 557 PetscCall(KSPGetType(ksp, &ksp_type)); 558 PetscCall(KSPGetConvergedReason(ksp, &reason)); 559 PetscCall(KSPGetIterationNumber(ksp, &its)); 560 PetscCall(KSPGetResidualNorm(ksp, &rnorm)); 561 562 if (!test_mode || reason < 0 || rnorm > 1e-8) { 563 PetscCall(PetscPrintf(comm, 564 "Swarm-to-Mesh Projection KSP Solve:\n" 565 " KSP type: %s\n" 566 " KSP convergence: %s\n" 567 " Total KSP iterations: %" PetscInt_FMT "\n" 568 " Final rnorm: %e\n", 569 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm)); 570 } 571 } 572 573 // Optional viewing 574 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_view")); 575 576 // Cleanup 577 PetscCall(VecDestroy(&B_mesh)); 578 PetscCall(MatDestroy(&M)); 579 PetscCall(KSPDestroy(&ksp)); 580 PetscFunctionReturn(PETSC_SUCCESS); 581 } 582 583 // ------------------------------------------------------------------------------------------------ 584 // BP setup 585 // ------------------------------------------------------------------------------------------------ 586 PetscErrorCode SetupProblemSwarm(DM dm_swarm, Ceed ceed, BPData bp_data, CeedData data, PetscBool setup_rhs, Vec rhs, Vec target) { 587 DM dm_mesh, dm_coord; 588 CeedElemRestriction elem_restr_u_mesh, elem_restr_x_mesh, elem_restr_x_points, elem_restr_u_points, elem_restr_q_data_points; 589 CeedBasis basis_u, basis_x; 590 CeedVector x_coord, x_ref_points, q_data_points; 591 PetscInt X_loc_len, dim; 592 CeedInt num_comp, q_data_size = bp_data.q_data_size; 593 CeedScalar R = 1; // radius of the sphere 594 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube 595 Vec X_loc; 596 PetscMemType X_mem_type; 597 598 PetscFunctionBeginUser; 599 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 600 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord)); 601 602 // Get coordinates 603 PetscCall(DMGetCoordinatesLocal(dm_mesh, &X_loc)); 604 PetscCall(VecGetLocalSize(X_loc, &X_loc_len)); 605 CeedVectorCreate(ceed, X_loc_len, &x_coord); 606 PetscCall(VecReadP2C(X_loc, &X_mem_type, x_coord)); 607 608 // Background mesh objects 609 PetscCall(CreateBasisFromPlex(ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &basis_u)); 610 PetscCall(CreateBasisFromPlex(ceed, dm_coord, NULL, 0, 0, 0, bp_data, &basis_x)); 611 PetscCall(CreateRestrictionFromPlex(ceed, dm_mesh, 0, NULL, 0, &elem_restr_u_mesh)); 612 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, NULL, 0, &elem_restr_x_mesh)); 613 614 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &data->x_ceed, NULL); 615 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &data->y_ceed, NULL); 616 617 // Swarm objects 618 { 619 const PetscInt *cell_points; 620 CeedInt *offsets; 621 IS is_points; 622 Vec X_ref; 623 CeedInt num_elem; 624 625 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref)); 626 PetscCall(DMGetDimension(dm_mesh, &dim)); 627 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem); 628 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp); 629 630 PetscCall(ISGetIndices(is_points, &cell_points)); 631 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2; 632 PetscCall(PetscCalloc1(num_elem + 1 + num_points, &offsets)); 633 634 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1; 635 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1]; 636 PetscCall(ISRestoreIndices(is_points, &cell_points)); 637 638 // -- Points restrictions 639 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 640 &elem_restr_u_points); 641 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 642 &elem_restr_x_points); 643 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, q_data_size, num_points * q_data_size, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 644 &elem_restr_q_data_points); 645 646 // -- Points vectors 647 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL); 648 649 // -- Ref coordinates 650 { 651 PetscMemType X_mem_type; 652 const PetscScalar *x; 653 654 CeedVectorCreate(ceed, num_points * dim, &x_ref_points); 655 656 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type)); 657 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x); 658 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x)); 659 } 660 661 // Create Q data 662 { 663 CeedQFunction qf_setup; 664 CeedOperator op_setup; 665 666 // Setup geometric scaling 667 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup); 668 CeedQFunctionAddInput(qf_setup, "x", dim, CEED_EVAL_INTERP); 669 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); 670 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 671 CeedQFunctionAddOutput(qf_setup, "qdata", q_data_size, CEED_EVAL_NONE); 672 673 CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 674 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 675 CeedOperatorSetField(op_setup, "dx", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 676 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 677 CeedOperatorSetField(op_setup, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 678 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points); 679 680 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE); 681 682 // Cleanup 683 CeedQFunctionDestroy(&qf_setup); 684 CeedOperatorDestroy(&op_setup); 685 } 686 687 // Cleanup 688 PetscCall(ISDestroy(&is_points)); 689 PetscCall(PetscFree(offsets)); 690 PetscCall(VecDestroy(&X_ref)); 691 } 692 693 // Set up PDE operator 694 695 CeedQFunction qf_apply; 696 CeedOperator op_apply; 697 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? dim : 1; 698 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? dim : 1; 699 700 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply); 701 CeedQFunctionAddInput(qf_apply, "u", num_comp * in_scale, bp_data.in_mode); 702 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 703 CeedQFunctionAddOutput(qf_apply, "v", num_comp * out_scale, bp_data.out_mode); 704 705 // Create the mass or diff operator 706 CeedOperatorCreateAtPoints(ceed, qf_apply, NULL, NULL, &op_apply); 707 CeedOperatorSetField(op_apply, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 708 CeedOperatorSetField(op_apply, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 709 CeedOperatorSetField(op_apply, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 710 CeedOperatorAtPointsSetPoints(op_apply, elem_restr_x_points, x_ref_points); 711 712 // Set up RHS if needed 713 if (setup_rhs) { 714 CeedQFunction qf_setup_rhs; 715 CeedOperator op_setup_rhs; 716 Vec rhs_loc; 717 CeedVector rhs_ceed, target_ceed; 718 PetscMemType rhs_mem_type, target_mem_type; 719 720 // Create RHS vector 721 PetscCall(DMCreateLocalVector(dm_mesh, &rhs_loc)); 722 723 CeedElemRestrictionCreateVector(elem_restr_u_points, &target_ceed, NULL); 724 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &rhs_ceed, NULL); 725 726 // Create the q-function that sets up the RHS and true solution 727 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 728 CeedQFunctionAddInput(qf_setup_rhs, "x", dim, CEED_EVAL_INTERP); 729 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 730 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp, CEED_EVAL_NONE); 731 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp, CEED_EVAL_INTERP); 732 733 // Create the operator that builds the RHS and true solution 734 CeedOperatorCreateAtPoints(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 735 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 736 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 737 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_points, CEED_BASIS_NONE, target_ceed); 738 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 739 CeedOperatorAtPointsSetPoints(op_setup_rhs, elem_restr_x_points, x_ref_points); 740 741 // Set up the libCEED context 742 CeedQFunctionContext ctx_rhs_setup; 743 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 744 CeedScalar rhs_setup_data[2] = {R, l}; 745 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 746 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 747 CeedQFunctionContextDestroy(&ctx_rhs_setup); 748 749 // Setup RHS and target 750 PetscCall(VecP2C(rhs_loc, &rhs_mem_type, rhs_ceed)); 751 PetscCall(VecP2C(target, &target_mem_type, target_ceed)); 752 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 753 PetscCall(VecC2P(rhs_ceed, rhs_mem_type, rhs_loc)); 754 PetscCall(VecC2P(target_ceed, target_mem_type, target)); 755 756 // Local-to-global 757 PetscCall(VecZeroEntries(rhs)); 758 PetscCall(DMLocalToGlobal(dm_mesh, rhs_loc, ADD_VALUES, rhs)); 759 760 PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view")); 761 762 // Cleanup 763 PetscCall(DMRestoreLocalVector(dm_mesh, &rhs_loc)); 764 CeedVectorDestroy(&rhs_ceed); 765 CeedVectorDestroy(&target_ceed); 766 CeedQFunctionDestroy(&qf_setup_rhs); 767 CeedOperatorDestroy(&op_setup_rhs); 768 } 769 770 // Save libCEED data 771 data->basis_x = basis_x; 772 data->basis_u = basis_u; 773 data->elem_restr_x = elem_restr_x_mesh; 774 data->elem_restr_u = elem_restr_u_mesh; 775 data->elem_restr_u_i = elem_restr_u_points; 776 data->elem_restr_qd_i = elem_restr_q_data_points; 777 data->qf_apply = qf_apply; 778 data->op_apply = op_apply; 779 data->q_data = q_data_points; 780 data->q_data_size = q_data_size; 781 782 // Cleanup 783 PetscCall(VecReadC2P(x_coord, X_mem_type, X_loc)); 784 CeedVectorDestroy(&x_coord); 785 CeedVectorDestroy(&x_ref_points); 786 CeedElemRestrictionDestroy(&elem_restr_x_points); 787 788 PetscFunctionReturn(PETSC_SUCCESS); 789 } 790