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(PetscFree(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 IS is_points; 621 Vec X_ref; 622 CeedInt num_elem; 623 624 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref)); 625 PetscCall(DMGetDimension(dm_mesh, &dim)); 626 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem); 627 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp); 628 629 PetscCall(ISGetIndices(is_points, &cell_points)); 630 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2; 631 CeedInt offsets[num_elem + 1 + num_points]; 632 633 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1; 634 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1]; 635 PetscCall(ISRestoreIndices(is_points, &cell_points)); 636 637 // -- Points restrictions 638 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 639 &elem_restr_u_points); 640 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 641 &elem_restr_x_points); 642 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, q_data_size, num_points * q_data_size, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 643 &elem_restr_q_data_points); 644 645 // -- Points vectors 646 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL); 647 648 // -- Ref coordinates 649 { 650 PetscMemType X_mem_type; 651 const PetscScalar *x; 652 653 CeedVectorCreate(ceed, num_points * dim, &x_ref_points); 654 655 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type)); 656 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x); 657 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x)); 658 } 659 660 // Create Q data 661 { 662 CeedQFunction qf_setup; 663 CeedOperator op_setup; 664 665 // Setup geometric scaling 666 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup); 667 CeedQFunctionAddInput(qf_setup, "x", dim, CEED_EVAL_INTERP); 668 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD); 669 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT); 670 CeedQFunctionAddOutput(qf_setup, "qdata", q_data_size, CEED_EVAL_NONE); 671 672 CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup); 673 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 674 CeedOperatorSetField(op_setup, "dx", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 675 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 676 CeedOperatorSetField(op_setup, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 677 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points); 678 679 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE); 680 681 // Cleanup 682 CeedQFunctionDestroy(&qf_setup); 683 CeedOperatorDestroy(&op_setup); 684 } 685 686 // Cleanup 687 PetscCall(ISDestroy(&is_points)); 688 PetscCall(VecDestroy(&X_ref)); 689 } 690 691 // Set up PDE operator 692 693 CeedQFunction qf_apply; 694 CeedOperator op_apply; 695 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? dim : 1; 696 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? dim : 1; 697 698 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply); 699 CeedQFunctionAddInput(qf_apply, "u", num_comp * in_scale, bp_data.in_mode); 700 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 701 CeedQFunctionAddOutput(qf_apply, "v", num_comp * out_scale, bp_data.out_mode); 702 703 // Create the mass or diff operator 704 CeedOperatorCreateAtPoints(ceed, qf_apply, NULL, NULL, &op_apply); 705 CeedOperatorSetField(op_apply, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 706 CeedOperatorSetField(op_apply, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 707 CeedOperatorSetField(op_apply, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 708 CeedOperatorAtPointsSetPoints(op_apply, elem_restr_x_points, x_ref_points); 709 710 // Set up RHS if needed 711 if (setup_rhs) { 712 CeedQFunction qf_setup_rhs; 713 CeedOperator op_setup_rhs; 714 Vec rhs_loc; 715 CeedVector rhs_ceed, target_ceed; 716 PetscMemType rhs_mem_type, target_mem_type; 717 718 // Create RHS vector 719 PetscCall(DMCreateLocalVector(dm_mesh, &rhs_loc)); 720 721 CeedElemRestrictionCreateVector(elem_restr_u_points, &target_ceed, NULL); 722 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &rhs_ceed, NULL); 723 724 // Create the q-function that sets up the RHS and true solution 725 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 726 CeedQFunctionAddInput(qf_setup_rhs, "x", dim, CEED_EVAL_INTERP); 727 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 728 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp, CEED_EVAL_NONE); 729 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp, CEED_EVAL_INTERP); 730 731 // Create the operator that builds the RHS and true solution 732 CeedOperatorCreateAtPoints(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 733 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE); 734 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points); 735 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_points, CEED_BASIS_NONE, target_ceed); 736 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE); 737 CeedOperatorAtPointsSetPoints(op_setup_rhs, elem_restr_x_points, x_ref_points); 738 739 // Set up the libCEED context 740 CeedQFunctionContext ctx_rhs_setup; 741 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 742 CeedScalar rhs_setup_data[2] = {R, l}; 743 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 744 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 745 CeedQFunctionContextDestroy(&ctx_rhs_setup); 746 747 // Setup RHS and target 748 PetscCall(VecP2C(rhs_loc, &rhs_mem_type, rhs_ceed)); 749 PetscCall(VecP2C(target, &target_mem_type, target_ceed)); 750 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 751 PetscCall(VecC2P(rhs_ceed, rhs_mem_type, rhs_loc)); 752 PetscCall(VecC2P(target_ceed, target_mem_type, target)); 753 754 // Local-to-global 755 PetscCall(VecZeroEntries(rhs)); 756 PetscCall(DMLocalToGlobal(dm_mesh, rhs_loc, ADD_VALUES, rhs)); 757 758 PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view")); 759 760 // Cleanup 761 PetscCall(DMRestoreLocalVector(dm_mesh, &rhs_loc)); 762 CeedVectorDestroy(&rhs_ceed); 763 CeedVectorDestroy(&target_ceed); 764 CeedQFunctionDestroy(&qf_setup_rhs); 765 CeedOperatorDestroy(&op_setup_rhs); 766 } 767 768 // Save libCEED data 769 data->basis_x = basis_x; 770 data->basis_u = basis_u; 771 data->elem_restr_x = elem_restr_x_mesh; 772 data->elem_restr_u = elem_restr_u_mesh; 773 data->elem_restr_u_i = elem_restr_u_points; 774 data->elem_restr_qd_i = elem_restr_q_data_points; 775 data->qf_apply = qf_apply; 776 data->op_apply = op_apply; 777 data->q_data = q_data_points; 778 data->q_data_size = q_data_size; 779 780 // Cleanup 781 PetscCall(VecReadC2P(x_coord, X_mem_type, X_loc)); 782 CeedVectorDestroy(&x_coord); 783 CeedVectorDestroy(&x_ref_points); 784 CeedElemRestrictionDestroy(&elem_restr_x_points); 785 786 PetscFunctionReturn(PETSC_SUCCESS); 787 } 788