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 // libCEED + PETSc DMSwarm Example 9 // 10 // This example demonstrates a simple usage of libCEED with DMSwarm. 11 // This example combines elements of PETSc src/impls/dm/swam/tutorials/ex1.c and src/impls/dm/swarm/tests/ex6.c 12 // 13 // Build with: 14 // 15 // make dmswarm [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 16 // 17 // Sample runs: 18 // 19 // ./dmswarm -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower -1.0,-1.0,-1.0 -dm_plex_simplex 0 -num_comp 2 -swarm gauss 20 // 21 //TESTARGS(name="Uniform swarm, CG projection") -ceed {ceed_resource} -test -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower -1.0,-1.0,-1.0 -dm_plex_simplex 0 -dm_plex_hash_location true -num_comp 2 -swarm uniform -solution_order 3 -points_per_cell 125 22 //TESTARGS(name="Gauss swarm, lumped projection") -ceed {ceed_resource} -test -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_box_lower -1.0,-1.0,-1.0 -dm_plex_simplex 0 -dm_plex_hash_location true -num_comp 2 -swarm gauss -ksp_type preonly -pc_type jacobi -pc_jacobi_type rowsum -tolerance 9e-2 23 24 /// @file 25 /// libCEED example using PETSc with DMSwarm 26 const char help[] = "libCEED example using PETSc with DMSwarm\n"; 27 28 #include <ceed.h> 29 #include <math.h> 30 #include <petscdmplex.h> 31 #include <petscdmswarm.h> 32 #include <petscds.h> 33 #include <petscfe.h> 34 #include <petscksp.h> 35 #include <petsc/private/petscfeimpl.h> /* For interpolation */ 36 37 #include "include/petscutils.h" 38 #include "include/petscversion.h" 39 #include "include/swarmutils.h" 40 41 const char DMSwarmPICField_u[] = "u"; 42 43 // Target functions 44 typedef enum { TARGET_TANH = 0, TARGET_POLYNOMIAL = 1, TARGET_SPHERE = 2 } TargetType; 45 static const char *const target_types[] = {"tanh", "polynomial", "sphere", "TargetType", "TARGET", 0}; 46 47 typedef PetscErrorCode (*TargetFunc)(PetscInt dim, const PetscScalar x[]); 48 typedef PetscErrorCode (*TargetFuncProj)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx); 49 50 PetscScalar EvalU_Tanh(PetscInt dim, const PetscScalar x[]); 51 PetscScalar EvalU_Poly(PetscInt dim, const PetscScalar x[]); 52 PetscScalar EvalU_Sphere(PetscInt dim, const PetscScalar x[]); 53 PetscErrorCode EvalU_Tanh_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 54 PetscErrorCode EvalU_Poly_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 55 PetscErrorCode EvalU_Sphere_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 56 57 // Swarm to mesh and mesh to swarm 58 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Petsc(DM dm_swarm, const char *field, Vec U_mesh); 59 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Ceed(DM dm_swarm, const char *field, Vec U_mesh); 60 PetscErrorCode DMSwarmCheckSwarmValues(DM dm_swarm, const char *field, PetscScalar tolerance, TargetFuncProj TrueSolution); 61 62 // ------------------------------------------------------------------------------------------------ 63 // main driver 64 // ------------------------------------------------------------------------------------------------ 65 int main(int argc, char **argv) { 66 MPI_Comm comm; 67 char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 68 PetscBool test_mode = PETSC_FALSE, view_petsc_swarm = PETSC_FALSE, view_ceed_swarm = PETSC_FALSE; 69 PetscInt dim = 3, num_comp = 1, num_points = 1728, num_points_per_cell = 64, mesh_order = 1, solution_order = 2, q_extra = 3; 70 PetscScalar tolerance = 1E-3; 71 DM dm_mesh, dm_swarm; 72 Vec U_mesh; 73 DMSwarmCeedContext swarm_ceed_context; 74 PointSwarmType point_swarm_type = SWARM_UNIFORM; 75 TargetType target_type = TARGET_TANH; 76 TargetFuncProj target_function_proj = EvalU_Tanh_proj; 77 78 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 79 comm = PETSC_COMM_WORLD; 80 81 // Read command line options 82 PetscOptionsBegin(comm, NULL, "libCEED example using PETSc with DMSwarm", NULL); 83 84 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 85 PetscCall( 86 PetscOptionsBool("-u_petsc_swarm_view", "View XDMF of swarm values interpolated by PETSc", NULL, view_petsc_swarm, &view_petsc_swarm, NULL)); 87 PetscCall( 88 PetscOptionsBool("-u_ceed_swarm_view", "View XDMF of swarm values interpolated by libCEED", NULL, view_ceed_swarm, &view_ceed_swarm, NULL)); 89 PetscCall(PetscOptionsEnum("-target", "Target field function", NULL, target_types, (PetscEnum)target_type, (PetscEnum *)&target_type, NULL)); 90 PetscCall(PetscOptionsInt("-solution_order", "Order of mesh solution space", NULL, solution_order, &solution_order, NULL)); 91 PetscCall(PetscOptionsInt("-mesh_order", "Order of mesh coordinate space", NULL, mesh_order, &mesh_order, NULL)); 92 PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, q_extra, &q_extra, NULL)); 93 PetscCall(PetscOptionsInt("-num_comp", "Number of components in solution", NULL, num_comp, &num_comp, NULL)); 94 PetscCall(PetscOptionsEnum("-swarm", "Swarm points distribution", NULL, point_swarm_types, (PetscEnum)point_swarm_type, 95 (PetscEnum *)&point_swarm_type, NULL)); 96 { 97 PetscBool user_set_num_points_per_cell = PETSC_FALSE; 98 PetscInt dim = 3, num_cells_total = 1; 99 PetscInt num_cells[] = {1, 1, 1}; 100 101 PetscCall(PetscOptionsInt("-points_per_cell", "Total number of swarm points in each cell", NULL, num_points_per_cell, &num_points_per_cell, 102 &user_set_num_points_per_cell)); 103 PetscCall(PetscOptionsInt("-dm_plex_dim", "Background mesh dimension", NULL, dim, &dim, NULL)); 104 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of cells", NULL, num_cells, &dim, NULL)); 105 106 num_cells_total = num_cells[0] * num_cells[1] * num_cells[2]; 107 PetscCheck(!user_set_num_points_per_cell || point_swarm_type != SWARM_SINUSOIDAL, comm, PETSC_ERR_USER, 108 "Cannot specify points per cell with sinusoidal points locations"); 109 if (!user_set_num_points_per_cell) { 110 PetscCall(PetscOptionsInt("-points", "Total number of swarm points", NULL, num_points, &num_points, NULL)); 111 num_points_per_cell = PetscCeilInt(num_points, num_cells_total); 112 } 113 if (point_swarm_type != SWARM_SINUSOIDAL) { 114 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)); 115 116 num_points_per_cell = 1; 117 for (PetscInt i = 0; i < dim; i++) num_points_per_cell *= num_points_per_cell_1d; 118 } 119 num_points = num_points_per_cell * num_cells_total; 120 } 121 PetscCall(PetscOptionsScalar("-tolerance", "Tolerance for swarm point values and projection relative L2 error", NULL, tolerance, &tolerance, NULL)); 122 PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, ceed_resource, ceed_resource, sizeof(ceed_resource), NULL)); 123 124 PetscOptionsEnd(); 125 126 // Create background mesh 127 { 128 PetscCall(DMCreate(comm, &dm_mesh)); 129 PetscCall(DMSetType(dm_mesh, DMPLEX)); 130 PetscCall(DMSetFromOptions(dm_mesh)); 131 132 // -- Check for tensor product mesh 133 { 134 PetscBool is_simplex; 135 136 PetscCall(DMPlexIsSimplex(dm_mesh, &is_simplex)); 137 PetscCheck(!is_simplex, comm, PETSC_ERR_USER, "Only tensor-product background meshes supported"); 138 } 139 140 // -- Mesh FE space 141 PetscCall(DMGetDimension(dm_mesh, &dim)); 142 { 143 PetscFE fe; 144 145 PetscCall(DMGetDimension(dm_mesh, &dim)); 146 PetscCall(PetscFECreateLagrange(comm, dim, num_comp, PETSC_FALSE, solution_order, solution_order + q_extra, &fe)); 147 PetscCall(DMAddField(dm_mesh, NULL, (PetscObject)fe)); 148 PetscCall(PetscFEDestroy(&fe)); 149 } 150 PetscCall(DMCreateDS(dm_mesh)); 151 152 // -- Coordinate FE space 153 { 154 PetscFE fe_coord; 155 156 PetscCall(PetscFECreateLagrange(comm, dim, dim, PETSC_FALSE, mesh_order, solution_order + q_extra, &fe_coord)); 157 PetscCall(DMSetCoordinateDisc(dm_mesh, fe_coord, PETSC_TRUE)); 158 PetscCall(PetscFEDestroy(&fe_coord)); 159 } 160 161 // -- Set tensor permutation 162 { 163 DM dm_coord; 164 165 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord)); 166 PetscCall(DMPlexSetClosurePermutationTensor(dm_mesh, PETSC_DETERMINE, NULL)); 167 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); 168 } 169 170 // -- Final background mesh 171 PetscCall(PetscObjectSetName((PetscObject)dm_mesh, "Background Mesh")); 172 PetscCall(DMViewFromOptions(dm_mesh, NULL, "-dm_mesh_view")); 173 } 174 175 // Create particle swarm 176 { 177 PetscCall(DMCreate(comm, &dm_swarm)); 178 PetscCall(DMSetType(dm_swarm, DMSWARM)); 179 PetscCall(DMSetDimension(dm_swarm, dim)); 180 PetscCall(DMSwarmSetType(dm_swarm, DMSWARM_PIC)); 181 PetscCall(DMSwarmSetCellDM(dm_swarm, dm_mesh)); 182 183 // -- Swarm field 184 PetscCall(DMSwarmRegisterPetscDatatypeField(dm_swarm, DMSwarmPICField_u, num_comp, PETSC_SCALAR)); 185 PetscCall(DMSwarmFinalizeFieldRegister(dm_swarm)); 186 PetscCall(DMSwarmSetLocalSizes(dm_swarm, num_points, 0)); 187 PetscCall(DMSetFromOptions(dm_swarm)); 188 189 // -- Set swarm point locations 190 PetscCall(DMSwarmInitalizePointLocations(dm_swarm, point_swarm_type, num_points, num_points_per_cell)); 191 192 // -- Final particle swarm 193 PetscCall(PetscObjectSetName((PetscObject)dm_swarm, "Particle Swarm")); 194 PetscCall(DMViewFromOptions(dm_swarm, NULL, "-dm_swarm_view")); 195 } 196 197 // Set field values on background mesh 198 PetscCall(DMCreateGlobalVector(dm_mesh, &U_mesh)); 199 switch (target_type) { 200 case TARGET_TANH: 201 target_function_proj = EvalU_Tanh_proj; 202 break; 203 case TARGET_POLYNOMIAL: 204 target_function_proj = EvalU_Poly_proj; 205 break; 206 case TARGET_SPHERE: 207 target_function_proj = EvalU_Sphere_proj; 208 break; 209 } 210 { 211 TargetFuncProj mesh_solution[1] = {target_function_proj}; 212 213 PetscCall(DMProjectFunction(dm_mesh, 0.0, mesh_solution, NULL, INSERT_VALUES, U_mesh)); 214 } 215 216 // Visualize background mesh 217 PetscCall(PetscObjectSetName((PetscObject)U_mesh, "U on Background Mesh")); 218 PetscCall(VecViewFromOptions(U_mesh, NULL, "-u_mesh_view")); 219 220 // libCEED objects for swarm and background mesh 221 PetscCall(DMSwarmCeedContextCreate(dm_swarm, ceed_resource, &swarm_ceed_context)); 222 223 // Interpolate from mesh to points via PETSc 224 { 225 PetscCall(DMSwarmInterpolateFromCellToSwarm_Petsc(dm_swarm, DMSwarmPICField_u, U_mesh)); 226 if (view_petsc_swarm) PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm_petsc.xmf")); 227 PetscCall(DMSwarmCheckSwarmValues(dm_swarm, DMSwarmPICField_u, tolerance, target_function_proj)); 228 } 229 230 // Interpolate from mesh to points via libCEED 231 { 232 PetscCall(DMSwarmInterpolateFromCellToSwarm_Ceed(dm_swarm, DMSwarmPICField_u, U_mesh)); 233 if (view_ceed_swarm) PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm_ceed.xmf")); 234 PetscCall(DMSwarmCheckSwarmValues(dm_swarm, DMSwarmPICField_u, tolerance, target_function_proj)); 235 } 236 237 // Project from points to mesh via libCEED 238 { 239 Vec U_projected; 240 241 PetscCall(VecDuplicate(U_mesh, &U_projected)); 242 PetscCall(DMSwarmProjectFromSwarmToCells(dm_swarm, DMSwarmPICField_u, U_projected)); 243 244 PetscCall(PetscObjectSetName((PetscObject)U_projected, "U projected to Background Mesh")); 245 PetscCall(VecViewFromOptions(U_projected, NULL, "-u_projected_view")); 246 PetscCall(VecAXPY(U_projected, -1.0, U_mesh)); 247 PetscCall(PetscObjectSetName((PetscObject)U_projected, "U Projection Error")); 248 PetscCall(VecViewFromOptions(U_projected, NULL, "-u_error_view")); 249 250 // -- Check error 251 { 252 PetscScalar error, norm_u_mesh; 253 254 PetscCall(VecNorm(U_projected, NORM_2, &error)); 255 PetscCall(VecNorm(U_mesh, NORM_2, &norm_u_mesh)); 256 PetscCheck(error / norm_u_mesh < tolerance, comm, PETSC_ERR_USER, "Projection error too high: %e\n", error / norm_u_mesh); 257 if (!test_mode) PetscCall(PetscPrintf(comm, " Projection error: %e\n", error / norm_u_mesh)); 258 } 259 260 PetscCall(VecDestroy(&U_projected)); 261 } 262 // Cleanup 263 PetscCall(DMSwarmCeedContextDestroy(&swarm_ceed_context)); 264 PetscCall(DMDestroy(&dm_swarm)); 265 PetscCall(DMDestroy(&dm_mesh)); 266 PetscCall(VecDestroy(&U_mesh)); 267 return PetscFinalize(); 268 } 269 270 // ------------------------------------------------------------------------------------------------ 271 // Solution functions 272 // ------------------------------------------------------------------------------------------------ 273 PetscScalar EvalU_Poly(PetscInt dim, const PetscScalar x[]) { 274 PetscScalar result = 0.0; 275 const PetscScalar p[5] = {3, 1, 4, 1, 5}; 276 277 for (PetscInt d = 0; d < dim; d++) { 278 PetscScalar result_1d = 1.0; 279 280 for (PetscInt i = 4; i >= 0; i--) result_1d = result_1d * x[d] + p[i]; 281 result += result_1d; 282 } 283 return result * 1E-3; 284 } 285 286 PetscScalar EvalU_Tanh(PetscInt dim, const PetscScalar x[]) { 287 PetscScalar result = 1.0, center = 0.1; 288 289 for (PetscInt d = 0; d < dim; d++) { 290 result *= tanh(x[d] - center); 291 center += 0.1; 292 } 293 return result; 294 } 295 296 PetscScalar EvalU_Sphere(PetscInt dim, const PetscScalar x[]) { 297 PetscScalar distance = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); 298 299 return distance < 1.0 ? 1.0 : 0.1; 300 } 301 302 PetscErrorCode EvalU_Poly_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 303 PetscFunctionBeginUser; 304 305 const PetscScalar f_x = EvalU_Poly(dim, x); 306 307 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 308 PetscFunctionReturn(PETSC_SUCCESS); 309 } 310 311 PetscErrorCode EvalU_Tanh_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 312 PetscFunctionBeginUser; 313 314 const PetscScalar f_x = EvalU_Tanh(dim, x); 315 316 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 317 PetscFunctionReturn(PETSC_SUCCESS); 318 } 319 320 PetscErrorCode EvalU_Sphere_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 321 PetscFunctionBeginUser; 322 323 const PetscScalar f_x = EvalU_Sphere(dim, x); 324 325 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 326 PetscFunctionReturn(PETSC_SUCCESS); 327 } 328 329 // ------------------------------------------------------------------------------------------------ 330 // Projection via PETSc 331 // ------------------------------------------------------------------------------------------------ 332 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Petsc(DM dm_swarm, const char *field, Vec U_mesh) { 333 PetscInt dim, num_comp, cell_start, cell_end; 334 PetscScalar *u_points; 335 const PetscScalar *coords_points; 336 const PetscReal v0_ref[3] = {-1.0, -1.0, -1.0}; 337 DM dm_mesh; 338 PetscSection section_u_mesh_loc; 339 PetscDS ds; 340 PetscFE fe; 341 PetscFEGeom fe_geometry; 342 PetscQuadrature quadrature; 343 Vec U_loc; 344 345 PetscFunctionBeginUser; 346 // Get mesh DM 347 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 348 PetscCall(DMGetDimension(dm_mesh, &dim)); 349 { 350 PetscSection section_u_mesh_loc_closure_permutation; 351 352 PetscCall(DMGetLocalSection(dm_mesh, §ion_u_mesh_loc_closure_permutation)); 353 PetscCall(PetscSectionClone(section_u_mesh_loc_closure_permutation, §ion_u_mesh_loc)); 354 PetscCall(PetscSectionResetClosurePermutation(section_u_mesh_loc)); 355 } 356 357 // Get local mesh values 358 PetscCall(DMGetLocalVector(dm_mesh, &U_loc)); 359 PetscCall(VecZeroEntries(U_loc)); 360 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_loc)); 361 362 // Get local swarm data 363 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 364 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 365 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 366 PetscCall(DMSwarmGetField(dm_swarm, field, &num_comp, NULL, (void **)&u_points)); 367 368 // Interpolate values to each swarm point, one element in the background mesh at a time 369 PetscCall(DMGetDS(dm_mesh, &ds)); 370 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe)); 371 for (PetscInt cell = cell_start; cell < cell_end; cell++) { 372 PetscTabulation tabulation; 373 PetscScalar *u_cell = NULL, *coords_points_cell_true, *coords_points_cell_ref; 374 PetscReal v[dim], J[dim * dim], invJ[dim * dim], detJ; 375 PetscInt *points_cell; 376 PetscInt num_points_in_cell; 377 378 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_cell)); 379 PetscCall(DMGetWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_true)); 380 PetscCall(DMGetWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_ref)); 381 // -- Reference coordinates for swarm points in background mesh element 382 for (PetscInt p = 0; p < num_points_in_cell; p++) { 383 for (PetscInt d = 0; d < dim; d++) coords_points_cell_true[p * dim + d] = coords_points[points_cell[p] * dim + d]; 384 } 385 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ)); 386 for (PetscInt p = 0; p < num_points_in_cell; p++) { 387 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_cell_true[p * dim], &coords_points_cell_ref[p * dim]); 388 } 389 // -- Interpolate values from current element in background mesh to swarm points 390 PetscCall(PetscFECreateTabulation(fe, 1, num_points_in_cell, coords_points_cell_ref, 1, &tabulation)); 391 PetscCall(DMPlexVecGetClosure(dm_mesh, section_u_mesh_loc, U_loc, cell, NULL, &u_cell)); 392 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 393 PetscCall(PetscFECreateCellGeometry(fe, quadrature, &fe_geometry)); 394 for (PetscInt p = 0; p < num_points_in_cell; p++) { 395 PetscCall(PetscFEInterpolateAtPoints_Static(fe, tabulation, u_cell, &fe_geometry, p, &u_points[points_cell[p] * num_comp])); 396 } 397 398 // -- Cleanup 399 PetscCall(PetscFEDestroyCellGeometry(fe, &fe_geometry)); 400 PetscCall(DMPlexVecRestoreClosure(dm_mesh, section_u_mesh_loc, U_loc, cell, NULL, &u_cell)); 401 PetscCall(DMRestoreWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_true)); 402 PetscCall(DMRestoreWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_ref)); 403 PetscCall(PetscTabulationDestroy(&tabulation)); 404 PetscCall(PetscFree(points_cell)); 405 } 406 407 // Cleanup 408 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 409 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&u_points)); 410 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 411 PetscCall(DMRestoreLocalVector(dm_mesh, &U_loc)); 412 PetscCall(PetscSectionDestroy(§ion_u_mesh_loc)); 413 PetscFunctionReturn(PETSC_SUCCESS); 414 } 415 416 // ------------------------------------------------------------------------------------------------ 417 // Projection via libCEED 418 // ------------------------------------------------------------------------------------------------ 419 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Ceed(DM dm_swarm, const char *field, Vec U_mesh) { 420 PetscMemType U_mem_type; 421 DM dm_mesh; 422 Vec U_mesh_loc; 423 DMSwarmCeedContext swarm_ceed_context; 424 425 PetscFunctionBeginUser; 426 // Get mesh DM 427 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 428 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 429 430 // Get mesh values 431 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc)); 432 PetscCall(VecZeroEntries(U_mesh_loc)); 433 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc)); 434 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh)); 435 436 // Get swarm access 437 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 438 PetscCall(DMSwarmPICFieldP2C(dm_swarm, field, swarm_ceed_context->u_points)); 439 440 // Interpolate field from mesh to swarm points 441 CeedOperatorApply(swarm_ceed_context->op_mesh_to_points, swarm_ceed_context->u_mesh, swarm_ceed_context->u_points, CEED_REQUEST_IMMEDIATE); 442 443 // Cleanup 444 PetscCall(DMSwarmPICFieldC2P(dm_swarm, field, swarm_ceed_context->u_points)); 445 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 446 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh, U_mem_type, U_mesh_loc)); 447 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc)); 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 // ------------------------------------------------------------------------------------------------ 452 // Error checking utility 453 // ------------------------------------------------------------------------------------------------ 454 PetscErrorCode DMSwarmCheckSwarmValues(DM dm_swarm, const char *field, PetscScalar tolerance, TargetFuncProj TrueSolution) { 455 PetscBool within_tolerance = PETSC_TRUE; 456 PetscInt dim, num_comp, cell_start, cell_end; 457 const PetscScalar *u_points, *coords_points; 458 DM dm_mesh; 459 460 PetscFunctionBeginUser; 461 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 462 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 463 PetscCall(DMGetDimension(dm_mesh, &dim)); 464 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 465 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 466 PetscCall(DMSwarmGetField(dm_swarm, field, &num_comp, NULL, (void **)&u_points)); 467 468 // Interpolate values to each swarm point, one element in the background mesh at a time 469 for (PetscInt cell = cell_start; cell < cell_end; cell++) { 470 PetscInt *points; 471 PetscInt num_points_in_cell; 472 473 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points)); 474 // -- Reference coordinates for swarm points in background mesh element 475 for (PetscInt p = 0; p < num_points_in_cell; p++) { 476 PetscScalar x[dim], u_true[num_comp]; 477 478 for (PetscInt d = 0; d < dim; d++) x[d] = coords_points[points[p] * dim + d]; 479 PetscCall(TrueSolution(dim, 0.0, x, num_comp, u_true, NULL)); 480 for (PetscInt i = 0; i < num_comp; i++) { 481 if (PetscAbs(u_points[points[p] * num_comp + i] - u_true[i]) > tolerance) { 482 within_tolerance = PETSC_FALSE; 483 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm_swarm), 484 "Incorrect interpolated value, cell %" PetscInt_FMT " point %" PetscInt_FMT " component %" PetscInt_FMT 485 ", found %f expected %f\n", 486 cell, p, i, u_points[points[p] * num_comp + i], u_true[i])); 487 } 488 } 489 } 490 491 // -- Cleanup 492 PetscCall(PetscFree(points)); 493 } 494 495 // Cleanup 496 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 497 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&u_points)); 498 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 499 PetscCheck(within_tolerance, PetscObjectComm((PetscObject)dm_swarm), PETSC_ERR_USER, "Interpolation to swarm points not within tolerance"); 500 PetscFunctionReturn(PETSC_SUCCESS); 501 } 502