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