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