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 40 const char DMSwarmPICField_u[] = "u"; 41 42 // libCEED context data 43 typedef struct DMSwarmCeedContext_ *DMSwarmCeedContext; 44 struct DMSwarmCeedContext_ { 45 Ceed ceed; 46 CeedVector u_mesh_loc, v_mesh_loc, u_mesh_elem, u_points_loc, u_points_elem, x_ref_points_loc, x_ref_points_elem; 47 CeedElemRestriction restriction_u_mesh, restriction_x_points, restriction_u_points; 48 CeedBasis basis_u; 49 }; 50 51 PetscErrorCode DMSwarmCeedContextCreate(DM dm_swarm, const char *ceed_resource, DMSwarmCeedContext *ctx); 52 PetscErrorCode DMSwarmCeedContextDestroy(DMSwarmCeedContext *ctx); 53 54 // Target functions 55 typedef enum { TARGET_TANH = 0, TARGET_POLYNOMIAL = 1, TARGET_SPHERE = 2 } TargetType; 56 static const char *const target_types[] = {"tanh", "polynomial", "sphere", "TargetType", "TARGET", 0}; 57 58 typedef PetscErrorCode (*TargetFunc)(PetscInt dim, const PetscScalar x[]); 59 typedef PetscErrorCode (*TargetFuncProj)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx); 60 61 PetscScalar EvalU_Tanh(PetscInt dim, const PetscScalar x[]); 62 PetscScalar EvalU_Poly(PetscInt dim, const PetscScalar x[]); 63 PetscScalar EvalU_Sphere(PetscInt dim, const PetscScalar x[]); 64 PetscErrorCode EvalU_Tanh_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 65 PetscErrorCode EvalU_Poly_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 66 PetscErrorCode EvalU_Sphere_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx); 67 68 // Swarm point distribution 69 typedef enum { SWARM_GAUSS = 0, SWARM_UNIFORM = 1, SWARM_CELL_RANDOM = 2, SWARM_SINUSOIDAL = 3 } PointSwarmType; 70 static const char *const point_swarm_types[] = {"gauss", "uniform", "cell_random", "sinusoidal", "PointSwarmType", "SWARM", 0}; 71 72 // Swarm helper function 73 PetscErrorCode DMSwarmInitalizePointLocations(DM dm_swarm, PointSwarmType point_swarm_type, PetscInt num_points, PetscInt num_points_per_cell); 74 75 // Memory utilities 76 PetscErrorCode VecP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed); 77 PetscErrorCode VecC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc); 78 PetscErrorCode VecReadP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed); 79 PetscErrorCode VecReadC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc); 80 PetscErrorCode DMSwarmPICFieldP2C(DM dm_swarm, const char *field, CeedVector x_ceed); 81 PetscErrorCode DMSwarmPICFieldC2P(DM dm_swarm, const char *field, CeedVector x_ceed); 82 83 // Swarm to mesh and mesh to swarm 84 PetscErrorCode DMSwarmCreateReferenceCoordinates(DM dm_swarm, IS *is_points, Vec *ref_coords); 85 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Petsc(DM dm_swarm, const char *field, Vec U_mesh); 86 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Ceed(DM dm_swarm, const char *field, Vec U_mesh); 87 PetscErrorCode DMSwarmCheckSwarmValues(DM dm_swarm, const char *field, PetscScalar tolerance, TargetFuncProj TrueSolution); 88 89 PetscErrorCode DMSwarmCreateProjectionRHS(DM dm_swarm, const char *field, Vec B_mesh); 90 PetscErrorCode MatMult_SwarmMass(Mat A, Vec U_mesh, Vec V_mesh); 91 92 // ------------------------------------------------------------------------------------------------ 93 // main driver 94 // ------------------------------------------------------------------------------------------------ 95 int main(int argc, char **argv) { 96 MPI_Comm comm; 97 char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; 98 PetscBool test_mode = PETSC_FALSE, view_petsc_swarm = PETSC_FALSE, view_ceed_swarm = PETSC_FALSE; 99 PetscInt dim = 3, num_comp = 1, num_points = 1728, num_points_per_cell = 64, mesh_order = 1, solution_order = 2, q_extra = 3; 100 PetscScalar tolerance = 1E-3; 101 DM dm_mesh, dm_swarm; 102 Vec U_mesh; 103 DMSwarmCeedContext swarm_ceed_context; 104 PointSwarmType point_swarm_type = SWARM_UNIFORM; 105 TargetType target_type = TARGET_TANH; 106 TargetFuncProj target_function_proj = EvalU_Tanh_proj; 107 108 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 109 comm = PETSC_COMM_WORLD; 110 111 // Read command line options 112 PetscOptionsBegin(comm, NULL, "libCEED example using PETSc with DMSwarm", NULL); 113 114 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 115 PetscCall( 116 PetscOptionsBool("-u_petsc_swarm_view", "View XDMF of swarm values interpolated by PETSc", NULL, view_petsc_swarm, &view_petsc_swarm, NULL)); 117 PetscCall( 118 PetscOptionsBool("-u_ceed_swarm_view", "View XDMF of swarm values interpolated by libCEED", NULL, view_ceed_swarm, &view_ceed_swarm, NULL)); 119 PetscCall(PetscOptionsEnum("-target", "Target field function", NULL, target_types, (PetscEnum)target_type, (PetscEnum *)&target_type, NULL)); 120 PetscCall(PetscOptionsInt("-solution_order", "Order of mesh solution space", NULL, solution_order, &solution_order, NULL)); 121 PetscCall(PetscOptionsInt("-mesh_order", "Order of mesh coordinate space", NULL, mesh_order, &mesh_order, NULL)); 122 PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, q_extra, &q_extra, NULL)); 123 PetscCall(PetscOptionsInt("-num_comp", "Number of components in solution", NULL, num_comp, &num_comp, NULL)); 124 PetscCall(PetscOptionsEnum("-swarm", "Swarm points distribution", NULL, point_swarm_types, (PetscEnum)point_swarm_type, 125 (PetscEnum *)&point_swarm_type, NULL)); 126 { 127 PetscBool user_set_num_points_per_cell = PETSC_FALSE; 128 PetscInt dim = 3, num_cells_total = 1; 129 PetscInt num_cells[] = {1, 1, 1}; 130 131 PetscCall(PetscOptionsInt("-points_per_cell", "Total number of swarm points in each cell", NULL, num_points_per_cell, &num_points_per_cell, 132 &user_set_num_points_per_cell)); 133 PetscCall(PetscOptionsInt("-dm_plex_dim", "Background mesh dimension", NULL, dim, &dim, NULL)); 134 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of cells", NULL, num_cells, &dim, NULL)); 135 136 num_cells_total = num_cells[0] * num_cells[1] * num_cells[2]; 137 PetscCheck(!user_set_num_points_per_cell || point_swarm_type != SWARM_SINUSOIDAL, comm, PETSC_ERR_USER, 138 "Cannot specify points per cell with sinusoidal points locations"); 139 if (!user_set_num_points_per_cell) { 140 PetscCall(PetscOptionsInt("-points", "Total number of swarm points", NULL, num_points, &num_points, NULL)); 141 num_points_per_cell = PetscCeilInt(num_points, num_cells_total); 142 } 143 if (point_swarm_type != SWARM_SINUSOIDAL) { 144 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)); 145 146 num_points_per_cell = 1; 147 for (PetscInt i = 0; i < dim; i++) num_points_per_cell *= num_points_per_cell_1d; 148 } 149 num_points = num_points_per_cell * num_cells_total; 150 } 151 PetscCall(PetscOptionsScalar("-tolerance", "Tolerance for swarm point values and projection relative L2 error", NULL, tolerance, &tolerance, NULL)); 152 PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, ceed_resource, ceed_resource, sizeof(ceed_resource), NULL)); 153 154 PetscOptionsEnd(); 155 156 // Create background mesh 157 { 158 PetscCall(DMCreate(comm, &dm_mesh)); 159 PetscCall(DMSetType(dm_mesh, DMPLEX)); 160 PetscCall(DMSetFromOptions(dm_mesh)); 161 162 // -- Check for tensor product mesh 163 { 164 PetscBool is_simplex; 165 166 PetscCall(DMPlexIsSimplex(dm_mesh, &is_simplex)); 167 PetscCheck(!is_simplex, comm, PETSC_ERR_USER, "Only tensor-product background meshes supported"); 168 } 169 170 // -- Mesh FE space 171 PetscCall(DMGetDimension(dm_mesh, &dim)); 172 { 173 PetscFE fe; 174 175 PetscCall(DMGetDimension(dm_mesh, &dim)); 176 PetscCall(PetscFECreateLagrange(comm, dim, num_comp, PETSC_FALSE, solution_order, solution_order + q_extra, &fe)); 177 PetscCall(DMAddField(dm_mesh, NULL, (PetscObject)fe)); 178 PetscCall(PetscFEDestroy(&fe)); 179 } 180 PetscCall(DMCreateDS(dm_mesh)); 181 182 // -- Coordinate FE space 183 { 184 PetscFE fe_coord; 185 186 PetscCall(PetscFECreateLagrange(comm, dim, dim, PETSC_FALSE, mesh_order, solution_order + q_extra, &fe_coord)); 187 PetscCall(DMSetCoordinateDisc(dm_mesh, fe_coord, PETSC_TRUE)); 188 PetscCall(PetscFEDestroy(&fe_coord)); 189 } 190 191 // -- Set tensor permutation 192 { 193 DM dm_coord; 194 195 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord)); 196 PetscCall(DMPlexSetClosurePermutationTensor(dm_mesh, PETSC_DETERMINE, NULL)); 197 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL)); 198 } 199 200 // -- Final background mesh 201 PetscCall(PetscObjectSetName((PetscObject)dm_mesh, "Background Mesh")); 202 PetscCall(DMViewFromOptions(dm_mesh, NULL, "-dm_mesh_view")); 203 } 204 205 // Create particle swarm 206 { 207 PetscCall(DMCreate(comm, &dm_swarm)); 208 PetscCall(DMSetType(dm_swarm, DMSWARM)); 209 PetscCall(DMSetDimension(dm_swarm, dim)); 210 PetscCall(DMSwarmSetType(dm_swarm, DMSWARM_PIC)); 211 PetscCall(DMSwarmSetCellDM(dm_swarm, dm_mesh)); 212 213 // -- Swarm field 214 PetscCall(DMSwarmRegisterPetscDatatypeField(dm_swarm, DMSwarmPICField_u, num_comp, PETSC_SCALAR)); 215 PetscCall(DMSwarmFinalizeFieldRegister(dm_swarm)); 216 PetscCall(DMSwarmSetLocalSizes(dm_swarm, num_points, 0)); 217 PetscCall(DMSetFromOptions(dm_swarm)); 218 219 // -- Set swarm point locations 220 PetscCall(DMSwarmInitalizePointLocations(dm_swarm, point_swarm_type, num_points, num_points_per_cell)); 221 222 // -- Final particle swarm 223 PetscCall(PetscObjectSetName((PetscObject)dm_swarm, "Particle Swarm")); 224 PetscCall(DMViewFromOptions(dm_swarm, NULL, "-dm_swarm_view")); 225 } 226 227 // Set field values on background mesh 228 PetscCall(DMCreateGlobalVector(dm_mesh, &U_mesh)); 229 switch (target_type) { 230 case TARGET_TANH: 231 target_function_proj = EvalU_Tanh_proj; 232 break; 233 case TARGET_POLYNOMIAL: 234 target_function_proj = EvalU_Poly_proj; 235 break; 236 case TARGET_SPHERE: 237 target_function_proj = EvalU_Sphere_proj; 238 break; 239 } 240 { 241 TargetFuncProj mesh_solution[1] = {target_function_proj}; 242 243 PetscCall(DMProjectFunction(dm_mesh, 0.0, mesh_solution, NULL, INSERT_VALUES, U_mesh)); 244 } 245 246 // Visualize background mesh 247 PetscCall(PetscObjectSetName((PetscObject)U_mesh, "U on Background Mesh")); 248 PetscCall(VecViewFromOptions(U_mesh, NULL, "-u_mesh_view")); 249 250 // libCEED objects for swarm and background mesh 251 PetscCall(DMSwarmCeedContextCreate(dm_swarm, ceed_resource, &swarm_ceed_context)); 252 253 // Interpolate from mesh to points via PETSc 254 { 255 PetscCall(DMSwarmInterpolateFromCellToSwarm_Petsc(dm_swarm, DMSwarmPICField_u, U_mesh)); 256 if (view_petsc_swarm) PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm_petsc.xmf")); 257 PetscCall(DMSwarmCheckSwarmValues(dm_swarm, DMSwarmPICField_u, tolerance, target_function_proj)); 258 } 259 260 // Interpolate from mesh to points via libCEED 261 { 262 PetscCall(DMSwarmInterpolateFromCellToSwarm_Ceed(dm_swarm, DMSwarmPICField_u, U_mesh)); 263 if (view_ceed_swarm) PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm_ceed.xmf")); 264 PetscCall(DMSwarmCheckSwarmValues(dm_swarm, DMSwarmPICField_u, tolerance, target_function_proj)); 265 } 266 267 // Project from points to mesh via libCEED 268 { 269 Vec B_mesh, U_projected; 270 Mat M; 271 KSP ksp; 272 273 PetscCall(VecDuplicate(U_mesh, &B_mesh)); 274 PetscCall(VecDuplicate(U_mesh, &U_projected)); 275 276 // -- Setup "mass matrix" 277 { 278 PetscInt l_size, g_size; 279 280 PetscCall(VecGetLocalSize(U_mesh, &l_size)); 281 PetscCall(VecGetSize(U_mesh, &g_size)); 282 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, swarm_ceed_context, &M)); 283 PetscCall(MatSetDM(M, dm_mesh)); 284 PetscCall(MatShellSetOperation(M, MATOP_MULT, (void (*)(void))MatMult_SwarmMass)); 285 } 286 287 // -- Setup KSP 288 { 289 PC pc; 290 291 PetscCall(KSPCreate(comm, &ksp)); 292 PetscCall(KSPGetPC(ksp, &pc)); 293 PetscCall(PCSetType(pc, PCJACOBI)); 294 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 295 PetscCall(KSPSetType(ksp, KSPCG)); 296 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 297 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 298 PetscCall(KSPSetOperators(ksp, M, M)); 299 PetscCall(KSPSetFromOptions(ksp)); 300 PetscCall(PetscObjectSetName((PetscObject)ksp, "Swarm-to-Mesh Projection")); 301 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_projection_view")); 302 } 303 304 // -- Setup RHS 305 PetscCall(DMSwarmCreateProjectionRHS(dm_swarm, DMSwarmPICField_u, B_mesh)); 306 307 // -- Solve 308 PetscCall(VecZeroEntries(U_projected)); 309 PetscCall(KSPSolve(ksp, B_mesh, U_projected)); 310 311 // -- KSP summary 312 { 313 KSPType ksp_type; 314 KSPConvergedReason reason; 315 PetscReal rnorm; 316 PetscInt its; 317 PetscCall(KSPGetType(ksp, &ksp_type)); 318 PetscCall(KSPGetConvergedReason(ksp, &reason)); 319 PetscCall(KSPGetIterationNumber(ksp, &its)); 320 PetscCall(KSPGetResidualNorm(ksp, &rnorm)); 321 322 if (!test_mode || reason < 0 || rnorm > 1e-8) { 323 PetscCall(PetscPrintf(comm, 324 "Swarm-to-Mesh Projection KSP Solve:\n" 325 " KSP type: %s\n" 326 " KSP convergence: %s\n" 327 " Total KSP iterations: %" PetscInt_FMT "\n" 328 " Final rnorm: %e\n", 329 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm)); 330 } 331 } 332 333 // -- Check error 334 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_view")); 335 PetscCall(PetscObjectSetName((PetscObject)U_mesh, "U projected to Background Mesh")); 336 PetscCall(VecViewFromOptions(U_projected, NULL, "-u_projected_view")); 337 PetscCall(VecAXPY(U_projected, -1.0, U_mesh)); 338 PetscCall(PetscObjectSetName((PetscObject)U_projected, "U Projection Error")); 339 PetscCall(VecViewFromOptions(U_projected, NULL, "-u_error_view")); 340 { 341 PetscScalar error, norm_u_mesh; 342 343 PetscCall(VecNorm(U_projected, NORM_2, &error)); 344 PetscCall(VecNorm(U_mesh, NORM_2, &norm_u_mesh)); 345 PetscCheck(error / norm_u_mesh < tolerance, comm, PETSC_ERR_USER, "Projection error too high: %e\n", error / norm_u_mesh); 346 if (!test_mode) PetscCall(PetscPrintf(comm, " Projection error: %e\n", error / norm_u_mesh)); 347 } 348 349 // -- Cleanup 350 PetscCall(VecDestroy(&B_mesh)); 351 PetscCall(VecDestroy(&U_projected)); 352 PetscCall(MatDestroy(&M)); 353 PetscCall(KSPDestroy(&ksp)); 354 } 355 356 // Cleanup 357 PetscCall(DMSwarmCeedContextDestroy(&swarm_ceed_context)); 358 PetscCall(DMDestroy(&dm_swarm)); 359 PetscCall(DMDestroy(&dm_mesh)); 360 PetscCall(VecDestroy(&U_mesh)); 361 return PetscFinalize(); 362 } 363 364 // ------------------------------------------------------------------------------------------------ 365 // Context utilities 366 // ------------------------------------------------------------------------------------------------ 367 PetscErrorCode DMSwarmCeedContextCreate(DM dm_swarm, const char *ceed_resource, DMSwarmCeedContext *ctx) { 368 DM dm_mesh; 369 370 PetscFunctionBeginUser; 371 PetscCall(PetscNew(ctx)); 372 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 373 374 CeedInit(ceed_resource, &(*ctx)->ceed); 375 // Background mesh objects 376 { 377 CeedInt elem_size, num_comp; 378 BPData bp_data = {.q_mode = CEED_GAUSS}; 379 380 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &(*ctx)->basis_u)); 381 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_mesh, 0, NULL, 0, &(*ctx)->restriction_u_mesh)); 382 383 // -- U vector 384 CeedElemRestrictionCreateVector((*ctx)->restriction_u_mesh, &(*ctx)->u_mesh_loc, NULL); 385 CeedElemRestrictionCreateVector((*ctx)->restriction_u_mesh, &(*ctx)->v_mesh_loc, NULL); 386 CeedElemRestrictionGetElementSize((*ctx)->restriction_u_mesh, &elem_size); 387 CeedElemRestrictionGetNumComponents((*ctx)->restriction_u_mesh, &num_comp); 388 CeedVectorCreate((*ctx)->ceed, elem_size * num_comp, &(*ctx)->u_mesh_elem); 389 } 390 // Swarm objects 391 { 392 PetscInt dim; 393 const PetscInt *cell_points; 394 IS is_points; 395 Vec X_ref; 396 CeedInt num_elem, num_comp, max_points_in_cell; 397 398 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref)); 399 PetscCall(DMGetDimension(dm_mesh, &dim)); 400 CeedElemRestrictionGetNumElements((*ctx)->restriction_u_mesh, &num_elem); 401 CeedElemRestrictionGetNumComponents((*ctx)->restriction_u_mesh, &num_comp); 402 403 PetscCall(ISGetIndices(is_points, &cell_points)); 404 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2; 405 CeedInt offsets[num_elem + 1 + num_points]; 406 407 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1; 408 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1]; 409 PetscCall(ISRestoreIndices(is_points, &cell_points)); 410 411 // -- Points restrictions 412 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 413 &(*ctx)->restriction_u_points); 414 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 415 &(*ctx)->restriction_x_points); 416 417 // -- U vector 418 CeedElemRestrictionGetMaxPointsInElement((*ctx)->restriction_u_points, &max_points_in_cell); 419 CeedElemRestrictionCreateVector((*ctx)->restriction_u_points, &(*ctx)->u_points_loc, NULL); 420 CeedVectorCreate((*ctx)->ceed, max_points_in_cell * num_comp, &(*ctx)->u_points_elem); 421 422 // -- Ref coordinates 423 { 424 PetscMemType X_mem_type; 425 const PetscScalar *x; 426 427 CeedVectorCreate((*ctx)->ceed, num_points * dim, &(*ctx)->x_ref_points_loc); 428 CeedVectorCreate((*ctx)->ceed, max_points_in_cell * dim, &(*ctx)->x_ref_points_elem); 429 430 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type)); 431 CeedVectorSetArray((*ctx)->x_ref_points_loc, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x); 432 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x)); 433 } 434 435 // -- Cleanup 436 PetscCall(ISDestroy(&is_points)); 437 PetscCall(VecDestroy(&X_ref)); 438 } 439 440 PetscCall(DMSetApplicationContext(dm_mesh, (void *)(*ctx))); 441 PetscFunctionReturn(PETSC_SUCCESS); 442 } 443 444 PetscErrorCode DMSwarmCeedContextDestroy(DMSwarmCeedContext *ctx) { 445 PetscFunctionBeginUser; 446 CeedDestroy(&(*ctx)->ceed); 447 CeedVectorDestroy(&(*ctx)->u_mesh_loc); 448 CeedVectorDestroy(&(*ctx)->v_mesh_loc); 449 CeedVectorDestroy(&(*ctx)->u_mesh_elem); 450 CeedVectorDestroy(&(*ctx)->u_points_loc); 451 CeedVectorDestroy(&(*ctx)->u_points_elem); 452 CeedVectorDestroy(&(*ctx)->x_ref_points_loc); 453 CeedVectorDestroy(&(*ctx)->x_ref_points_elem); 454 CeedElemRestrictionDestroy(&(*ctx)->restriction_u_mesh); 455 CeedElemRestrictionDestroy(&(*ctx)->restriction_x_points); 456 CeedElemRestrictionDestroy(&(*ctx)->restriction_u_points); 457 CeedBasisDestroy(&(*ctx)->basis_u); 458 PetscCall(PetscFree(*ctx)); 459 PetscFunctionReturn(PETSC_SUCCESS); 460 } 461 462 // ------------------------------------------------------------------------------------------------ 463 // PETSc-libCEED memory space utilities 464 // ------------------------------------------------------------------------------------------------ 465 PetscErrorCode VecP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) { 466 PetscScalar *x; 467 468 PetscFunctionBeginUser; 469 PetscCall(VecGetArrayAndMemType(X_petsc, &x, mem_type)); 470 CeedVectorSetArray(x_ceed, MemTypeP2C(*mem_type), CEED_USE_POINTER, x); 471 PetscFunctionReturn(PETSC_SUCCESS); 472 } 473 474 PetscErrorCode VecC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) { 475 PetscScalar *x; 476 477 PetscFunctionBeginUser; 478 CeedVectorTakeArray(x_ceed, MemTypeP2C(mem_type), &x); 479 PetscCall(VecRestoreArrayAndMemType(X_petsc, &x)); 480 PetscFunctionReturn(PETSC_SUCCESS); 481 } 482 483 PetscErrorCode VecReadP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) { 484 PetscScalar *x; 485 486 PetscFunctionBeginUser; 487 PetscCall(VecGetArrayReadAndMemType(X_petsc, (const PetscScalar **)&x, mem_type)); 488 CeedVectorSetArray(x_ceed, MemTypeP2C(*mem_type), CEED_USE_POINTER, x); 489 PetscFunctionReturn(PETSC_SUCCESS); 490 } 491 492 PetscErrorCode VecReadC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) { 493 PetscScalar *x; 494 495 PetscFunctionBeginUser; 496 CeedVectorTakeArray(x_ceed, MemTypeP2C(mem_type), &x); 497 PetscCall(VecRestoreArrayReadAndMemType(X_petsc, (const PetscScalar **)&x)); 498 PetscFunctionReturn(PETSC_SUCCESS); 499 } 500 501 PetscErrorCode DMSwarmPICFieldP2C(DM dm_swarm, const char *field, CeedVector x_ceed) { 502 PetscScalar *x; 503 504 PetscFunctionBeginUser; 505 PetscCall(DMSwarmGetField(dm_swarm, field, NULL, NULL, (void **)&x)); 506 CeedVectorSetArray(x_ceed, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)x); 507 PetscFunctionReturn(PETSC_SUCCESS); 508 } 509 510 PetscErrorCode DMSwarmPICFieldC2P(DM dm_swarm, const char *field, CeedVector x_ceed) { 511 PetscScalar *x; 512 513 PetscFunctionBeginUser; 514 CeedVectorTakeArray(x_ceed, CEED_MEM_HOST, (CeedScalar **)&x); 515 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&x)); 516 PetscFunctionReturn(PETSC_SUCCESS); 517 } 518 519 // ------------------------------------------------------------------------------------------------ 520 // Solution functions 521 // ------------------------------------------------------------------------------------------------ 522 PetscScalar EvalU_Poly(PetscInt dim, const PetscScalar x[]) { 523 PetscScalar result = 0.0; 524 const PetscScalar p[5] = {3, 1, 4, 1, 5}; 525 526 for (PetscInt d = 0; d < dim; d++) { 527 PetscScalar result_1d = 1.0; 528 529 for (PetscInt i = 4; i >= 0; i--) result_1d = result_1d * x[d] + p[i]; 530 result += result_1d; 531 } 532 return result * 1E-3; 533 } 534 535 PetscScalar EvalU_Tanh(PetscInt dim, const PetscScalar x[]) { 536 PetscScalar result = 1.0, center = 0.1; 537 538 for (PetscInt d = 0; d < dim; d++) { 539 result *= tanh(x[d] - center); 540 center += 0.1; 541 } 542 return result; 543 } 544 545 PetscScalar EvalU_Sphere(PetscInt dim, const PetscScalar x[]) { 546 PetscScalar distance = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]); 547 548 return distance < 1.0 ? 1.0 : 0.1; 549 } 550 551 PetscErrorCode EvalU_Poly_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 552 PetscFunctionBeginUser; 553 554 const PetscScalar f_x = EvalU_Poly(dim, x); 555 556 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 557 PetscFunctionReturn(PETSC_SUCCESS); 558 } 559 560 PetscErrorCode EvalU_Tanh_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 561 PetscFunctionBeginUser; 562 563 const PetscScalar f_x = EvalU_Tanh(dim, x); 564 565 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 566 PetscFunctionReturn(PETSC_SUCCESS); 567 } 568 569 PetscErrorCode EvalU_Sphere_proj(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt num_comp, PetscScalar *u, void *ctx) { 570 PetscFunctionBeginUser; 571 572 const PetscScalar f_x = EvalU_Sphere(dim, x); 573 574 for (PetscInt c = 0; c < num_comp; c++) u[c] = (c + 1.0) * f_x; 575 PetscFunctionReturn(PETSC_SUCCESS); 576 } 577 578 // ------------------------------------------------------------------------------------------------ 579 // Swarm point location utility 580 // ------------------------------------------------------------------------------------------------ 581 PetscErrorCode DMSwarmInitalizePointLocations(DM dm_swarm, PointSwarmType point_swarm_type, PetscInt num_points, PetscInt num_points_per_cell) { 582 PetscFunctionBeginUser; 583 584 switch (point_swarm_type) { 585 case SWARM_GAUSS: 586 case SWARM_UNIFORM: { 587 // -- Set gauss or uniform point locations in each cell 588 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)), dim = 3; 589 PetscScalar point_coords[num_points_per_cell * 3]; 590 CeedScalar points_1d[num_points_per_cell_1d], weights_1d[num_points_per_cell_1d]; 591 592 if (point_swarm_type == SWARM_GAUSS) { 593 PetscCall(CeedGaussQuadrature(num_points_per_cell_1d, points_1d, weights_1d)); 594 } else { 595 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) points_1d[i] = 2.0 * (PetscReal)(i + 1) / (PetscReal)(num_points_per_cell_1d + 1) - 1; 596 } 597 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) { 598 for (PetscInt j = 0; j < num_points_per_cell_1d; j++) { 599 for (PetscInt k = 0; k < num_points_per_cell_1d; k++) { 600 PetscInt p = (i * num_points_per_cell_1d + j) * num_points_per_cell_1d + k; 601 602 point_coords[p * dim + 0] = points_1d[i]; 603 point_coords[p * dim + 1] = points_1d[j]; 604 point_coords[p * dim + 2] = points_1d[k]; 605 } 606 } 607 } 608 PetscCall(DMSwarmSetPointCoordinatesCellwise(dm_swarm, num_points_per_cell_1d * num_points_per_cell_1d * num_points_per_cell_1d, point_coords)); 609 } break; 610 case SWARM_CELL_RANDOM: { 611 // -- Set points randomly in each cell 612 PetscInt dim = 3, num_cells_total = 1, num_cells[] = {1, 1, 1}; 613 PetscScalar *point_coords; 614 PetscRandom rng; 615 616 PetscOptionsBegin(PetscObjectComm((PetscObject)dm_swarm), NULL, "libCEED example using PETSc with DMSwarm", NULL); 617 618 PetscCall(PetscOptionsInt("-dm_plex_dim", "Background mesh dimension", NULL, dim, &dim, NULL)); 619 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of cells", NULL, num_cells, &dim, NULL)); 620 621 PetscOptionsEnd(); 622 623 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm_swarm), &rng)); 624 625 num_cells_total = num_cells[0] * num_cells[1] * num_cells[2]; 626 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 627 for (PetscInt c = 0; c < num_cells_total; c++) { 628 PetscInt cell_index[3] = {c % num_cells[0], (c / num_cells[0]) % num_cells[1], (c / num_cells[0] / num_cells[1]) % num_cells[2]}; 629 630 for (PetscInt p = 0; p < num_points_per_cell; p++) { 631 PetscInt point_index = c * num_points_per_cell + p; 632 PetscScalar random_value; 633 634 for (PetscInt i = 0; i < dim; i++) { 635 PetscCall(PetscRandomGetValue(rng, &random_value)); 636 point_coords[point_index * dim + i] = -1.0 + cell_index[i] * 2.0 / (num_cells[i] + 1.0) + random_value; 637 } 638 } 639 } 640 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 641 PetscCall(PetscRandomDestroy(&rng)); 642 } break; 643 case SWARM_SINUSOIDAL: { 644 // -- Set points distributed per sinusoidal functions 645 PetscInt dim = 3; 646 PetscScalar *point_coords; 647 DM dm_mesh; 648 649 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 650 PetscCall(DMGetDimension(dm_mesh, &dim)); 651 652 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 653 for (PetscInt p = 0; p < num_points; p++) { 654 point_coords[p * dim + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 655 if (dim > 1) point_coords[p * dim + 1] = -PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 656 if (dim > 2) point_coords[p * dim + 2] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI); 657 } 658 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords)); 659 } break; 660 } 661 PetscCall(DMSwarmMigrate(dm_swarm, PETSC_TRUE)); 662 PetscFunctionReturn(PETSC_SUCCESS); 663 } 664 665 /*@C 666 DMSwarmCreateReferenceCoordinates - Compute the cell reference coordinates for local DMSwarm points. 667 668 Collective 669 670 Input Parameter: 671 . dm_swarm - the `DMSwarm` 672 673 Output Parameters: 674 + is_points - The IS object for indexing into points per cell 675 - X_points_ref - Vec holding the cell reference coordinates for local DMSwarm points 676 677 The index set contains ranges of indices for each local cell. This range contains the indices of every point in the cell. 678 679 ``` 680 total_num_cells 681 cell_0_start_index 682 cell_1_start_index 683 cell_2_start_index 684 ... 685 cell_n_start_index 686 cell_n_stop_index 687 cell_0_point_0 688 cell_0_point_0 689 ... 690 cell_n_point_m 691 ``` 692 693 Level: beginner 694 695 .seealso: `DMSwarm` 696 @*/ 697 PetscErrorCode DMSwarmCreateReferenceCoordinates(DM dm_swarm, IS *is_points, Vec *X_points_ref) { 698 PetscInt cell_start, cell_end, num_cells_local, dim, num_points_local, *cell_points, points_offset; 699 PetscScalar *coords_points_ref; 700 const PetscScalar *coords_points_true; 701 DM dm_mesh; 702 703 PetscFunctionBeginUser; 704 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 705 706 // Create vector to hold reference coordinates 707 { 708 Vec X_points_true; 709 710 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 711 PetscCall(VecDuplicate(X_points_true, X_points_ref)); 712 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true)); 713 } 714 715 // Allocate index set array 716 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 717 num_cells_local = cell_end - cell_start; 718 points_offset = num_cells_local + 2; 719 PetscCall(VecGetLocalSize(*X_points_ref, &num_points_local)); 720 PetscCall(DMGetDimension(dm_mesh, &dim)); 721 num_points_local /= dim; 722 PetscCall(PetscMalloc1(num_points_local + num_cells_local + 2, &cell_points)); 723 cell_points[0] = num_cells_local; 724 725 // Get reference coordinates for each swarm point wrt the elements in the background mesh 726 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 727 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 728 PetscCall(VecGetArray(*X_points_ref, &coords_points_ref)); 729 for (PetscInt cell = cell_start, num_points_processed = 0; cell < cell_end; cell++) { 730 PetscInt *points_in_cell, num_points_in_cell, local_cell = cell - cell_start; 731 PetscReal v[3], J[9], invJ[9], detJ, v0_ref[3] = {-1.0, -1.0, -1.0}; 732 733 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell)); 734 // -- Reference coordinates for swarm points in background mesh element 735 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ)); 736 cell_points[local_cell + 1] = num_points_processed + points_offset; 737 for (PetscInt p = 0; p < num_points_in_cell; p++) { 738 const CeedInt point = points_in_cell[p]; 739 740 cell_points[points_offset + (num_points_processed++)] = point; 741 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_true[point * dim], &coords_points_ref[point * dim]); 742 } 743 744 // -- Cleanup 745 PetscCall(PetscFree(points_in_cell)); 746 } 747 cell_points[points_offset - 1] = num_points_local + points_offset; 748 749 // Cleanup 750 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true)); 751 PetscCall(VecRestoreArray(*X_points_ref, &coords_points_ref)); 752 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 753 754 // Create index set 755 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_points_local + points_offset, cell_points, PETSC_OWN_POINTER, is_points)); 756 PetscFunctionReturn(PETSC_SUCCESS); 757 } 758 759 // ------------------------------------------------------------------------------------------------ 760 // Projection via PETSc 761 // ------------------------------------------------------------------------------------------------ 762 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Petsc(DM dm_swarm, const char *field, Vec U_mesh) { 763 PetscInt dim, num_comp, cell_start, cell_end; 764 PetscScalar *u_points; 765 const PetscScalar *coords_points; 766 const PetscReal v0_ref[3] = {-1.0, -1.0, -1.0}; 767 DM dm_mesh; 768 PetscSection section_u_mesh_loc; 769 PetscDS ds; 770 PetscFE fe; 771 PetscFEGeom fe_geometry; 772 PetscQuadrature quadrature; 773 Vec U_loc; 774 775 PetscFunctionBeginUser; 776 // Get mesh DM 777 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 778 PetscCall(DMGetDimension(dm_mesh, &dim)); 779 { 780 PetscSection section_u_mesh_loc_closure_permutation; 781 782 PetscCall(DMGetLocalSection(dm_mesh, §ion_u_mesh_loc_closure_permutation)); 783 PetscCall(PetscSectionClone(section_u_mesh_loc_closure_permutation, §ion_u_mesh_loc)); 784 PetscCall(PetscSectionResetClosurePermutation(section_u_mesh_loc)); 785 } 786 787 // Get local mesh values 788 PetscCall(DMGetLocalVector(dm_mesh, &U_loc)); 789 PetscCall(VecZeroEntries(U_loc)); 790 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_loc)); 791 792 // Get local swarm data 793 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 794 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 795 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 796 PetscCall(DMSwarmGetField(dm_swarm, field, &num_comp, NULL, (void **)&u_points)); 797 798 // Interpolate values to each swarm point, one element in the background mesh at a time 799 PetscCall(DMGetDS(dm_mesh, &ds)); 800 PetscCall(PetscDSGetDiscretization(ds, 0, (PetscObject *)&fe)); 801 for (PetscInt cell = cell_start; cell < cell_end; cell++) { 802 PetscTabulation tabulation; 803 PetscScalar *u_cell = NULL, *coords_points_cell_true, *coords_points_cell_ref; 804 PetscReal v[dim], J[dim * dim], invJ[dim * dim], detJ; 805 PetscInt *points_cell; 806 PetscInt num_points_in_cell; 807 808 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_cell)); 809 PetscCall(DMGetWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_true)); 810 PetscCall(DMGetWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_ref)); 811 // -- Reference coordinates for swarm points in background mesh element 812 for (PetscInt p = 0; p < num_points_in_cell; p++) { 813 for (PetscInt d = 0; d < dim; d++) coords_points_cell_true[p * dim + d] = coords_points[points_cell[p] * dim + d]; 814 } 815 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ)); 816 for (PetscInt p = 0; p < num_points_in_cell; p++) { 817 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_cell_true[p * dim], &coords_points_cell_ref[p * dim]); 818 } 819 // -- Interpolate values from current element in background mesh to swarm points 820 PetscCall(PetscFECreateTabulation(fe, 1, num_points_in_cell, coords_points_cell_ref, 1, &tabulation)); 821 PetscCall(DMPlexVecGetClosure(dm_mesh, section_u_mesh_loc, U_loc, cell, NULL, &u_cell)); 822 PetscCall(PetscFEGetQuadrature(fe, &quadrature)); 823 PetscCall(PetscFECreateCellGeometry(fe, quadrature, &fe_geometry)); 824 for (PetscInt p = 0; p < num_points_in_cell; p++) { 825 PetscCall(PetscFEInterpolateAtPoints_Static(fe, tabulation, u_cell, &fe_geometry, p, &u_points[points_cell[p] * num_comp])); 826 } 827 828 // -- Cleanup 829 PetscCall(PetscFEDestroyCellGeometry(fe, &fe_geometry)); 830 PetscCall(DMPlexVecRestoreClosure(dm_mesh, section_u_mesh_loc, U_loc, cell, NULL, &u_cell)); 831 PetscCall(DMRestoreWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_true)); 832 PetscCall(DMRestoreWorkArray(dm_mesh, num_points_in_cell * dim, MPIU_REAL, &coords_points_cell_ref)); 833 PetscCall(PetscTabulationDestroy(&tabulation)); 834 PetscCall(PetscFree(points_cell)); 835 } 836 837 // Cleanup 838 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 839 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&u_points)); 840 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 841 PetscCall(DMRestoreLocalVector(dm_mesh, &U_loc)); 842 PetscCall(PetscSectionDestroy(§ion_u_mesh_loc)); 843 PetscFunctionReturn(PETSC_SUCCESS); 844 } 845 846 // ------------------------------------------------------------------------------------------------ 847 // Projection via libCEED 848 // ------------------------------------------------------------------------------------------------ 849 PetscErrorCode DMSwarmInterpolateFromCellToSwarm_Ceed(DM dm_swarm, const char *field, Vec U_mesh) { 850 PetscInt num_elem; 851 PetscMemType U_mem_type; 852 DM dm_mesh; 853 Vec U_mesh_loc; 854 DMSwarmCeedContext swarm_ceed_context; 855 856 PetscFunctionBeginUser; 857 // Get mesh DM 858 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 859 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 860 861 // Get mesh values 862 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc)); 863 PetscCall(VecZeroEntries(U_mesh_loc)); 864 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc)); 865 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh_loc)); 866 867 // Get swarm access 868 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 869 PetscCall(DMSwarmPICFieldP2C(dm_swarm, field, swarm_ceed_context->u_points_loc)); 870 871 // Interpolate values to each swarm point, one element in the background mesh at a time 872 CeedElemRestrictionGetNumElements(swarm_ceed_context->restriction_u_mesh, &num_elem); 873 for (PetscInt e = 0; e < num_elem; e++) { 874 PetscInt num_points_in_elem; 875 876 CeedElemRestrictionGetNumPointsInElement(swarm_ceed_context->restriction_u_points, e, &num_points_in_elem); 877 878 // -- Reference coordinates for swarm points in background mesh element 879 CeedElemRestrictionApplyAtPointsInElement(swarm_ceed_context->restriction_x_points, e, CEED_NOTRANSPOSE, swarm_ceed_context->x_ref_points_loc, 880 swarm_ceed_context->x_ref_points_elem, CEED_REQUEST_IMMEDIATE); 881 882 // -- Interpolate values from current element in background mesh to swarm points 883 // Note: This will only work for CPU backends at this time, as only CPU backends support ApplyBlock and ApplyAtPoints 884 CeedElemRestrictionApplyBlock(swarm_ceed_context->restriction_u_mesh, e, CEED_NOTRANSPOSE, swarm_ceed_context->u_mesh_loc, 885 swarm_ceed_context->u_mesh_elem, CEED_REQUEST_IMMEDIATE); 886 CeedBasisApplyAtPoints(swarm_ceed_context->basis_u, num_points_in_elem, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, swarm_ceed_context->x_ref_points_elem, 887 swarm_ceed_context->u_mesh_elem, swarm_ceed_context->u_points_elem); 888 889 // -- Insert result back into local vector 890 CeedElemRestrictionApplyAtPointsInElement(swarm_ceed_context->restriction_u_points, e, CEED_TRANSPOSE, swarm_ceed_context->u_points_elem, 891 swarm_ceed_context->u_points_loc, CEED_REQUEST_IMMEDIATE); 892 } 893 894 // Cleanup 895 PetscCall(DMSwarmPICFieldC2P(dm_swarm, field, swarm_ceed_context->u_points_loc)); 896 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 897 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh_loc, U_mem_type, U_mesh_loc)); 898 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc)); 899 PetscFunctionReturn(PETSC_SUCCESS); 900 } 901 902 // ------------------------------------------------------------------------------------------------ 903 // Error checking utility 904 // ------------------------------------------------------------------------------------------------ 905 PetscErrorCode DMSwarmCheckSwarmValues(DM dm_swarm, const char *field, PetscScalar tolerance, TargetFuncProj TrueSolution) { 906 PetscBool within_tolerance = PETSC_TRUE; 907 PetscInt dim, num_comp, cell_start, cell_end; 908 const PetscScalar *u_points, *coords_points; 909 DM dm_mesh; 910 911 PetscFunctionBeginUser; 912 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 913 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 914 PetscCall(DMGetDimension(dm_mesh, &dim)); 915 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end)); 916 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 917 PetscCall(DMSwarmGetField(dm_swarm, field, &num_comp, NULL, (void **)&u_points)); 918 919 // Interpolate values to each swarm point, one element in the background mesh at a time 920 for (PetscInt cell = cell_start; cell < cell_end; cell++) { 921 PetscInt *points; 922 PetscInt num_points_in_cell; 923 924 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points)); 925 // -- Reference coordinates for swarm points in background mesh element 926 for (PetscInt p = 0; p < num_points_in_cell; p++) { 927 PetscScalar x[dim], u_true[num_comp]; 928 929 for (PetscInt d = 0; d < dim; d++) x[d] = coords_points[points[p] * dim + d]; 930 PetscCall(TrueSolution(dim, 0.0, x, num_comp, u_true, NULL)); 931 for (PetscInt i = 0; i < num_comp; i++) { 932 if (PetscAbs(u_points[points[p] * num_comp + i] - u_true[i]) > tolerance) { 933 within_tolerance = PETSC_FALSE; 934 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm_swarm), 935 "Incorrect interpolated value, cell %" PetscInt_FMT " point %" PetscInt_FMT " component %" PetscInt_FMT 936 ", found %f expected %f\n", 937 cell, p, i, u_points[points[p] * num_comp + i], u_true[i])); 938 } 939 } 940 } 941 942 // -- Cleanup 943 PetscCall(PetscFree(points)); 944 } 945 946 // Cleanup 947 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points)); 948 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&u_points)); 949 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 950 PetscCheck(within_tolerance, PetscObjectComm((PetscObject)dm_swarm), PETSC_ERR_USER, "Interpolation to swarm points not within tolerance"); 951 PetscFunctionReturn(PETSC_SUCCESS); 952 } 953 954 // ------------------------------------------------------------------------------------------------ 955 // RHS for Swarm to Mesh projection 956 // ------------------------------------------------------------------------------------------------ 957 PetscErrorCode DMSwarmCreateProjectionRHS(DM dm_swarm, const char *field, Vec B_mesh) { 958 PetscInt num_elem; 959 PetscMemType B_mem_type; 960 DM dm_mesh; 961 Vec B_mesh_loc; 962 DMSwarmCeedContext swarm_ceed_context; 963 964 PetscFunctionBeginUser; 965 // Get mesh DM 966 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh)); 967 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 968 969 // Get mesh values 970 PetscCall(DMGetLocalVector(dm_mesh, &B_mesh_loc)); 971 PetscCall(VecZeroEntries(B_mesh_loc)); 972 PetscCall(VecP2C(B_mesh_loc, &B_mem_type, swarm_ceed_context->v_mesh_loc)); 973 974 // Get swarm access 975 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 976 PetscCall(DMSwarmPICFieldP2C(dm_swarm, field, swarm_ceed_context->u_points_loc)); 977 978 // Interpolate values to each swarm point, one element in the background mesh at a time 979 CeedElemRestrictionGetNumElements(swarm_ceed_context->restriction_u_mesh, &num_elem); 980 for (PetscInt e = 0; e < num_elem; e++) { 981 PetscInt num_points_in_elem; 982 983 CeedElemRestrictionGetNumPointsInElement(swarm_ceed_context->restriction_u_points, e, &num_points_in_elem); 984 985 // -- Reference coordinates for swarm points in background mesh element 986 CeedElemRestrictionApplyAtPointsInElement(swarm_ceed_context->restriction_x_points, e, CEED_NOTRANSPOSE, swarm_ceed_context->x_ref_points_loc, 987 swarm_ceed_context->x_ref_points_elem, CEED_REQUEST_IMMEDIATE); 988 989 // -- Interpolate values from current element in background mesh to swarm points 990 // Note: This will only work for CPU backends at this time, as only CPU backends support ApplyBlock and ApplyAtPoints 991 CeedElemRestrictionApplyAtPointsInElement(swarm_ceed_context->restriction_u_points, e, CEED_NOTRANSPOSE, swarm_ceed_context->u_points_loc, 992 swarm_ceed_context->u_points_elem, CEED_REQUEST_IMMEDIATE); 993 CeedBasisApplyAtPoints(swarm_ceed_context->basis_u, num_points_in_elem, CEED_TRANSPOSE, CEED_EVAL_INTERP, swarm_ceed_context->x_ref_points_elem, 994 swarm_ceed_context->u_points_elem, swarm_ceed_context->u_mesh_elem); 995 996 // -- Insert result back into local vector 997 CeedElemRestrictionApplyBlock(swarm_ceed_context->restriction_u_mesh, e, CEED_TRANSPOSE, swarm_ceed_context->u_mesh_elem, 998 swarm_ceed_context->v_mesh_loc, CEED_REQUEST_IMMEDIATE); 999 } 1000 1001 // Restore PETSc Vecs and Local to Global 1002 PetscCall(VecC2P(swarm_ceed_context->v_mesh_loc, B_mem_type, B_mesh_loc)); 1003 PetscCall(VecZeroEntries(B_mesh)); 1004 PetscCall(DMLocalToGlobal(dm_mesh, B_mesh_loc, ADD_VALUES, B_mesh)); 1005 1006 // Cleanup 1007 PetscCall(DMSwarmPICFieldC2P(dm_swarm, field, swarm_ceed_context->u_points_loc)); 1008 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 1009 PetscCall(DMRestoreLocalVector(dm_mesh, &B_mesh_loc)); 1010 PetscFunctionReturn(PETSC_SUCCESS); 1011 } 1012 1013 // ------------------------------------------------------------------------------------------------ 1014 // Swarm "mass matrix" 1015 // ------------------------------------------------------------------------------------------------ 1016 PetscErrorCode MatMult_SwarmMass(Mat A, Vec U_mesh, Vec V_mesh) { 1017 PetscInt num_elem; 1018 PetscMemType U_mem_type, V_mem_type; 1019 DM dm_mesh; 1020 Vec U_mesh_loc, V_mesh_loc; 1021 DMSwarmCeedContext swarm_ceed_context; 1022 1023 PetscFunctionBeginUser; 1024 // Get mesh DM 1025 PetscCall(MatGetDM(A, &dm_mesh)); 1026 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context)); 1027 1028 // Global to Local and get PETSc Vec access 1029 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc)); 1030 PetscCall(VecZeroEntries(U_mesh_loc)); 1031 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc)); 1032 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh_loc)); 1033 PetscCall(DMGetLocalVector(dm_mesh, &V_mesh_loc)); 1034 PetscCall(VecZeroEntries(V_mesh_loc)); 1035 PetscCall(VecP2C(V_mesh_loc, &V_mem_type, swarm_ceed_context->v_mesh_loc)); 1036 1037 // Interpolate values to each swarm point, one element in the background mesh at a time 1038 CeedElemRestrictionGetNumElements(swarm_ceed_context->restriction_u_mesh, &num_elem); 1039 for (PetscInt e = 0; e < num_elem; e++) { 1040 PetscInt num_points_in_elem; 1041 1042 CeedElemRestrictionGetNumPointsInElement(swarm_ceed_context->restriction_u_points, e, &num_points_in_elem); 1043 1044 // -- Reference coordinates for swarm points in background mesh element 1045 CeedElemRestrictionApplyAtPointsInElement(swarm_ceed_context->restriction_x_points, e, CEED_NOTRANSPOSE, swarm_ceed_context->x_ref_points_loc, 1046 swarm_ceed_context->x_ref_points_elem, CEED_REQUEST_IMMEDIATE); 1047 1048 // -- Interpolate values from current element in background mesh to swarm points 1049 // Note: This will only work for CPU backends at this time, as only CPU backends support ApplyBlock and ApplyAtPoints 1050 CeedElemRestrictionApplyBlock(swarm_ceed_context->restriction_u_mesh, e, CEED_NOTRANSPOSE, swarm_ceed_context->u_mesh_loc, 1051 swarm_ceed_context->u_mesh_elem, CEED_REQUEST_IMMEDIATE); 1052 CeedBasisApplyAtPoints(swarm_ceed_context->basis_u, num_points_in_elem, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, swarm_ceed_context->x_ref_points_elem, 1053 swarm_ceed_context->u_mesh_elem, swarm_ceed_context->u_points_elem); 1054 1055 // -- Interpolate transpose back from swarm points to mesh 1056 CeedBasisApplyAtPoints(swarm_ceed_context->basis_u, num_points_in_elem, CEED_TRANSPOSE, CEED_EVAL_INTERP, swarm_ceed_context->x_ref_points_elem, 1057 swarm_ceed_context->u_points_elem, swarm_ceed_context->u_mesh_elem); 1058 CeedElemRestrictionApplyBlock(swarm_ceed_context->restriction_u_mesh, e, CEED_TRANSPOSE, swarm_ceed_context->u_mesh_elem, 1059 swarm_ceed_context->v_mesh_loc, CEED_REQUEST_IMMEDIATE); 1060 } 1061 1062 // Restore PETSc Vecs and Local to Global 1063 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh_loc, U_mem_type, U_mesh_loc)); 1064 PetscCall(VecC2P(swarm_ceed_context->v_mesh_loc, V_mem_type, V_mesh_loc)); 1065 PetscCall(VecZeroEntries(V_mesh)); 1066 PetscCall(DMLocalToGlobal(dm_mesh, V_mesh_loc, ADD_VALUES, V_mesh)); 1067 1068 // Cleanup 1069 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc)); 1070 PetscCall(DMRestoreLocalVector(dm_mesh, &V_mesh_loc)); 1071 PetscFunctionReturn(PETSC_SUCCESS); 1072 } 1073