1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 // libCEED + PETSc Example: CEED BPs 9 // 10 // This example demonstrates a simple usage of libCEED with PETSc to solve the CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps, on 11 // a particle swarm. 12 // 13 // The code uses higher level communication protocols in DMPlex and DMSwarm. 14 // 15 // Build with: 16 // 17 // make bpsswarm [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 18 // 19 // Sample runs: 20 // 21 // bpssphere -problem bp1 -degree 3 22 // bpssphere -problem bp2 -degree 3 23 // bpssphere -problem bp3 -degree 3 24 // 25 //TESTARGS(name="BP2") -ceed {ceed_resource} -test -problem bp2 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_simplex 0 -swarm uniform -points_per_cell 64 26 //TESTARGS(name="BP3") -ceed {ceed_resource} -test -problem bp3 -dm_plex_dim 3 -dm_plex_box_faces 4,4,4 -dm_plex_simplex 0 -swarm uniform -points_per_cell 64 -tolerance 3e-2 27 //TESTARGS(name="BP5") -ceed {ceed_resource} -test -problem bp5 -dm_plex_dim 3 -dm_plex_box_faces 3,3,3 -dm_plex_simplex 0 -swarm uniform -points_per_cell 64 28 29 /// @file 30 /// CEED BPs example using PETSc with DMPlex 31 /// See bpsraw.c for a "raw" implementation using a structured grid and bps.c for an implementation using an unstructured grid. 32 static const char help[] = "Solve CEED BPs on a particle swarm using DMPlex and DMSwarm in PETSc\n"; 33 const char DMSwarmPICField_u[] = "u"; 34 35 #include "bps.h" 36 37 #include <ceed.h> 38 #include <petscdmplex.h> 39 #include <petscksp.h> 40 #include <stdbool.h> 41 #include <string.h> 42 43 #include "include/bpsproblemdata.h" 44 #include "include/libceedsetup.h" 45 #include "include/matops.h" 46 #include "include/petscutils.h" 47 #include "include/petscversion.h" 48 #include "include/swarmutils.h" 49 50 int main(int argc, char **argv) { 51 MPI_Comm comm; 52 char ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self", filename[PETSC_MAX_PATH_LEN]; 53 double my_rt_start, my_rt, rt_min, rt_max; 54 PetscScalar tolerance; 55 PetscMPIInt comm_size; 56 PetscInt degree, q_extra, l_size, g_size, dim = 3, num_comp_u = 1, xl_size, num_points = 1728, num_points_per_cell = 64; 57 PetscBool test_mode, benchmark_mode, read_mesh, write_solution, write_true_solution_swarm; 58 PetscLogStage solve_stage; 59 Vec X, X_loc, rhs; 60 Mat mat_O; 61 KSP ksp; 62 DM dm_mesh, dm_swarm; 63 OperatorApplyContext op_apply_ctx, op_error_ctx; 64 Ceed ceed; 65 CeedData ceed_data; 66 CeedOperator op_error; 67 BPType bp_choice; 68 VecType vec_type; 69 PointSwarmType point_swarm_type = SWARM_GAUSS; 70 PetscMPIInt ranks_per_node; 71 char hostname[PETSC_MAX_PATH_LEN]; 72 73 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 74 comm = PETSC_COMM_WORLD; 75 PetscCall(MPI_Comm_size(comm, &comm_size)); 76 #if defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY) 77 { 78 MPI_Comm splitcomm; 79 PetscCall(MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &splitcomm)); 80 PetscCall(MPI_Comm_size(splitcomm, &ranks_per_node)); 81 PetscCall(MPI_Comm_free(&splitcomm)); 82 } 83 #else 84 ranks_per_node = -1; // Unknown 85 #endif 86 87 // Read command line options 88 PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); 89 bp_choice = CEED_BP1; 90 PetscCall(PetscOptionsEnum("-problem", "CEED benchmark problem to solve", NULL, bp_types, (PetscEnum)bp_choice, (PetscEnum *)&bp_choice, NULL)); 91 num_comp_u = bp_options[bp_choice].num_comp_u; 92 test_mode = PETSC_FALSE; 93 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL)); 94 benchmark_mode = PETSC_FALSE; 95 PetscCall(PetscOptionsBool("-benchmark", "Benchmarking mode (prints benchmark statistics)", NULL, benchmark_mode, &benchmark_mode, NULL)); 96 write_solution = PETSC_FALSE; 97 PetscCall(PetscOptionsBool("-write_solution", "Write solution for visualization", NULL, write_solution, &write_solution, NULL)); 98 write_true_solution_swarm = PETSC_FALSE; 99 PetscCall(PetscOptionsBool("-write_true_solution_swarm", "Write true solution at swarm points for visualization", NULL, write_true_solution_swarm, 100 &write_true_solution_swarm, NULL)); 101 degree = 2; 102 PetscCall(PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", NULL, degree, °ree, NULL)); 103 q_extra = bp_options[bp_choice].q_extra; 104 PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, q_extra, &q_extra, NULL)); 105 PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, ceed_resource, ceed_resource, sizeof(ceed_resource), NULL)); 106 PetscCall(PetscGetHostName(hostname, sizeof hostname)); 107 PetscCall(PetscOptionsString("-hostname", "Hostname for output", NULL, hostname, hostname, sizeof(hostname), NULL)); 108 read_mesh = PETSC_FALSE; 109 PetscCall(PetscOptionsString("-mesh", "Read mesh from file", NULL, filename, filename, sizeof(filename), &read_mesh)); 110 tolerance = 1e-2; 111 PetscCall(PetscOptionsScalar("-tolerance", "Tolerance for L2 error", NULL, tolerance, &tolerance, NULL)); 112 PetscCall(PetscOptionsEnum("-swarm", "Swarm points distribution", NULL, point_swarm_types, (PetscEnum)point_swarm_type, 113 (PetscEnum *)&point_swarm_type, NULL)); 114 { 115 PetscBool user_set_num_points_per_cell = PETSC_FALSE; 116 PetscInt num_cells_total = 1, tmp = dim; 117 PetscInt num_cells[] = {1, 1, 1}; 118 119 PetscCall(PetscOptionsInt("-points_per_cell", "Total number of swarm points in each cell", NULL, num_points_per_cell, &num_points_per_cell, 120 &user_set_num_points_per_cell)); 121 PetscCall(PetscOptionsInt("-dm_plex_dim", "Background mesh dimension", NULL, dim, &dim, NULL)); 122 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of cells", NULL, num_cells, &tmp, NULL)); 123 124 PetscCheck(tmp == dim, comm, PETSC_ERR_USER, "Number of values for -dm_plex_box_faces must match dimension"); 125 126 num_cells_total = num_cells[0] * num_cells[1] * num_cells[2]; 127 PetscCheck(!user_set_num_points_per_cell || point_swarm_type != SWARM_SINUSOIDAL, comm, PETSC_ERR_USER, 128 "Cannot specify points per cell with sinusoidal points locations"); 129 if (!user_set_num_points_per_cell) { 130 PetscCall(PetscOptionsInt("-points", "Total number of swarm points", NULL, num_points, &num_points, NULL)); 131 num_points_per_cell = PetscCeilInt(num_points, num_cells_total); 132 } 133 if (point_swarm_type != SWARM_SINUSOIDAL) { 134 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)); 135 136 num_points_per_cell = 1; 137 for (PetscInt i = 0; i < dim; i++) num_points_per_cell *= num_points_per_cell_1d; 138 } 139 num_points = num_points_per_cell * num_cells_total; 140 } 141 { 142 PetscBool flg; 143 PetscInt p = ranks_per_node; 144 PetscCall(PetscOptionsInt("-p", "Number of MPI ranks per node", NULL, p, &p, &flg)); 145 if (flg) ranks_per_node = p; 146 } 147 PetscOptionsEnd(); 148 149 // Setup DM 150 if (read_mesh) { 151 PetscCall(DMPlexCreateFromFile(comm, filename, NULL, PETSC_TRUE, &dm_mesh)); 152 } else { 153 PetscCall(DMCreate(comm, &dm_mesh)); 154 PetscCall(DMSetType(dm_mesh, DMPLEX)); 155 PetscCall(DMSetFromOptions(dm_mesh)); 156 157 // -- Check for tensor product mesh 158 { 159 PetscBool is_simplex; 160 161 PetscCall(DMPlexIsSimplex(dm_mesh, &is_simplex)); 162 PetscCheck(!is_simplex, comm, PETSC_ERR_USER, "Only tensor-product background meshes supported"); 163 } 164 } 165 PetscCall(DMGetDimension(dm_mesh, &dim)); 166 PetscCall(SetupDMByDegree(dm_mesh, degree, q_extra, num_comp_u, dim, bp_options[bp_choice].enforce_bc)); 167 168 // View mesh 169 PetscCall(DMSetOptionsPrefix(dm_mesh, "final_")); 170 PetscCall(DMViewFromOptions(dm_mesh, NULL, "-dm_view")); 171 172 // Create particle swarm 173 PetscCall(DMCreate(comm, &dm_swarm)); 174 PetscCall(DMSetType(dm_swarm, DMSWARM)); 175 PetscCall(DMSetDimension(dm_swarm, dim)); 176 PetscCall(DMSwarmSetType(dm_swarm, DMSWARM_PIC)); 177 PetscCall(DMSwarmSetCellDM(dm_swarm, dm_mesh)); 178 179 // -- Swarm field 180 PetscCall(DMSwarmRegisterPetscDatatypeField(dm_swarm, DMSwarmPICField_u, num_comp_u, PETSC_SCALAR)); 181 PetscCall(DMSwarmFinalizeFieldRegister(dm_swarm)); 182 { 183 PetscInt c_start, c_end, num_cells_local; 184 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &c_start, &c_end)); 185 num_cells_local = c_end - c_start; 186 PetscCall(DMSwarmSetLocalSizes(dm_swarm, num_cells_local * num_points_per_cell, 0)); 187 } 188 PetscCall(DMSetFromOptions(dm_swarm)); 189 190 // -- Set swarm point locations 191 PetscCall(DMSwarmInitalizePointLocations(dm_swarm, point_swarm_type, num_points, num_points_per_cell)); 192 PetscCall(DMSwarmVectorDefineField(dm_swarm, DMSwarmPICField_u)); 193 194 // -- Final particle swarm 195 PetscCall(PetscObjectSetName((PetscObject)dm_swarm, "Particle Swarm")); 196 PetscCall(DMViewFromOptions(dm_swarm, NULL, "-dm_swarm_view")); 197 198 // Create vectors 199 PetscCall(DMCreateGlobalVector(dm_mesh, &X)); 200 PetscCall(VecGetLocalSize(X, &l_size)); 201 PetscCall(VecGetSize(X, &g_size)); 202 PetscCall(DMCreateLocalVector(dm_mesh, &X_loc)); 203 PetscCall(VecGetSize(X_loc, &xl_size)); 204 PetscCall(VecDuplicate(X, &rhs)); 205 206 // Operator 207 PetscCall(PetscMalloc1(1, &op_apply_ctx)); 208 PetscCall(PetscMalloc1(1, &op_error_ctx)); 209 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, op_apply_ctx, &mat_O)); 210 PetscCall(MatSetDM(mat_O, dm_mesh)); 211 PetscCall(MatShellSetOperation(mat_O, MATOP_MULT, (void (*)(void))MatMult_Ceed)); 212 213 // Set up libCEED 214 CeedInit(ceed_resource, &ceed); 215 CeedMemType mem_type_backend; 216 CeedGetPreferredMemType(ceed, &mem_type_backend); 217 218 PetscCall(DMGetVecType(dm_mesh, &vec_type)); 219 if (!vec_type) { // Not yet set by user -dm_vec_type 220 switch (mem_type_backend) { 221 case CEED_MEM_HOST: 222 vec_type = VECSTANDARD; 223 break; 224 case CEED_MEM_DEVICE: { 225 const char *resolved; 226 CeedGetResource(ceed, &resolved); 227 if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; 228 else if (strstr(resolved, "/gpu/hip/occa")) vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 229 else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP; 230 else vec_type = VECSTANDARD; 231 } 232 } 233 PetscCall(DMSetVecType(dm_mesh, vec_type)); 234 } 235 236 // Print summary 237 if (!test_mode) { 238 PetscInt P = degree + 1, Q = P + q_extra; 239 240 const char *used_resource; 241 CeedGetResource(ceed, &used_resource); 242 243 VecType vec_type; 244 PetscCall(VecGetType(X, &vec_type)); 245 246 PetscInt c_start, c_end, num_cells_local; 247 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &c_start, &c_end)); 248 num_cells_local = c_end - c_start; 249 DMPolytopeType cell_type; 250 PetscCall(DMPlexGetCellType(dm_mesh, c_start, &cell_type)); 251 PetscMPIInt comm_size; 252 PetscCall(MPI_Comm_size(comm, &comm_size)); 253 254 PetscInt num_points_local, num_points_global; 255 PetscCall(DMSwarmGetLocalSize(dm_swarm, &num_points_local)); 256 PetscCall(DMSwarmGetSize(dm_swarm, &num_points_global)); 257 258 PetscCall(PetscPrintf(comm, 259 "\n-- CEED Benchmark Problem %" CeedInt_FMT " -- libCEED + PETSc --\n" 260 " MPI:\n" 261 " Hostname : %s\n" 262 " Total ranks : %d\n" 263 " Ranks per compute node : %d\n" 264 " PETSc:\n" 265 " PETSc Vec Type : %s\n" 266 " libCEED:\n" 267 " libCEED Backend : %s\n" 268 " libCEED Backend MemType : %s\n" 269 " Mesh:\n" 270 " Solution Order (P) : %" PetscInt_FMT "\n" 271 " Quadrature Order (Q) : %" PetscInt_FMT "\n" 272 " Additional quadrature points (q_extra) : %" PetscInt_FMT "\n" 273 " Global nodes : %" PetscInt_FMT "\n" 274 " Local Elements : %" PetscInt_FMT "\n" 275 " Owned nodes : %" PetscInt_FMT "\n" 276 " DoF per node : %" PetscInt_FMT "\n" 277 " Swarm:\n" 278 " Global points : %" PetscInt_FMT "\n" 279 " Local points : %" PetscInt_FMT "\n" 280 " Avg points per cell : %" PetscInt_FMT "\n" 281 " Point distribution : %s\n", 282 bp_choice + 1, hostname, comm_size, ranks_per_node, vec_type, used_resource, CeedMemTypes[mem_type_backend], P, Q, q_extra, 283 g_size / num_comp_u, num_cells_local, l_size / num_comp_u, num_comp_u, num_points_global, num_points_local, 284 num_cells_local > 0 ? num_points_local / num_cells_local : 0, point_swarm_types[point_swarm_type])); 285 } 286 287 // Setup libCEED's objects 288 Vec target; 289 290 PetscCall(DMCreateLocalVector(dm_swarm, &target)); 291 PetscCall(PetscMalloc1(1, &ceed_data)); 292 PetscCall(SetupProblemSwarm(dm_swarm, ceed, bp_options[bp_choice], ceed_data, true, rhs, target)); 293 PetscCall(SetupErrorOperator(dm_mesh, ceed, bp_options[bp_choice], dim, dim, num_comp_u, &op_error)); 294 295 // Set up apply operator context 296 PetscCall(SetupApplyOperatorCtx(comm, dm_mesh, ceed, ceed_data, X_loc, op_apply_ctx)); 297 298 // Setup solver 299 PetscCall(KSPCreate(comm, &ksp)); 300 { 301 PC pc; 302 PetscCall(KSPGetPC(ksp, &pc)); 303 if (bp_choice == CEED_BP1 || bp_choice == CEED_BP2) { 304 PetscCall(PCSetType(pc, PCJACOBI)); 305 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM)); 306 } else { 307 PetscCall(PCSetType(pc, PCNONE)); 308 MatNullSpace nullspace; 309 310 PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace)); 311 PetscCall(MatSetNullSpace(mat_O, nullspace)); 312 PetscCall(MatNullSpaceDestroy(&nullspace)); 313 } 314 PetscCall(KSPSetType(ksp, KSPCG)); 315 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 316 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 317 } 318 PetscCall(KSPSetFromOptions(ksp)); 319 PetscCall(KSPSetOperators(ksp, mat_O, mat_O)); 320 321 // First run, if benchmarking 322 if (benchmark_mode) { 323 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1)); 324 my_rt_start = MPI_Wtime(); 325 PetscCall(KSPSolve(ksp, rhs, X)); 326 my_rt = MPI_Wtime() - my_rt_start; 327 PetscCall(MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm)); 328 // Set maxits based on first iteration timing 329 if (my_rt > 0.02) { 330 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 5)); 331 } else { 332 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 20)); 333 } 334 } 335 336 // Timed solve 337 PetscCall(VecZeroEntries(X)); 338 PetscCall(PetscBarrier((PetscObject)ksp)); 339 340 // -- Performance logging 341 PetscCall(PetscLogStageRegister("Solve Stage", &solve_stage)); 342 PetscCall(PetscLogStagePush(solve_stage)); 343 344 // -- Solve 345 my_rt_start = MPI_Wtime(); 346 PetscCall(KSPSolve(ksp, rhs, X)); 347 my_rt = MPI_Wtime() - my_rt_start; 348 349 // -- Performance logging 350 PetscCall(PetscLogStagePop()); 351 352 // Output results 353 { 354 KSPType ksp_type; 355 KSPConvergedReason reason; 356 PetscReal rnorm; 357 PetscInt its; 358 PetscCall(KSPGetType(ksp, &ksp_type)); 359 PetscCall(KSPGetConvergedReason(ksp, &reason)); 360 PetscCall(KSPGetIterationNumber(ksp, &its)); 361 PetscCall(KSPGetResidualNorm(ksp, &rnorm)); 362 if (!test_mode || reason < 0 || rnorm > 1e-8) { 363 PetscCall(PetscPrintf(comm, 364 " KSP:\n" 365 " KSP Type : %s\n" 366 " KSP Convergence : %s\n" 367 " Total KSP Iterations : %" PetscInt_FMT "\n" 368 " Final rnorm : %e\n", 369 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm)); 370 } 371 if (!test_mode) { 372 PetscCall(PetscPrintf(comm, " Performance:\n")); 373 } 374 375 // View true solution at particles 376 if (write_true_solution_swarm) { 377 Vec u_swarm, u_swarm_old; 378 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 379 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm)); 380 PetscCall(VecDuplicate(u_swarm, &u_swarm_old)); 381 PetscCall(VecCopy(u_swarm, u_swarm_old)); 382 PetscCall(VecCopy(target, u_swarm)); 383 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm)); 384 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 385 386 PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm.xmf")); 387 388 PetscCall(DMSwarmSortGetAccess(dm_swarm)); 389 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm)); 390 PetscCall(VecCopy(u_swarm_old, u_swarm)); 391 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm)); 392 PetscCall(DMSwarmSortRestoreAccess(dm_swarm)); 393 PetscCall(VecDestroy(&u_swarm_old)); 394 } 395 396 // View solution at mesh points 397 PetscCall(VecViewFromOptions(X, NULL, "-solution_view")); 398 399 // Compute L2 Error 400 { 401 // Set up error operator context 402 PetscCall(SetupErrorOperatorCtx(comm, dm_mesh, ceed, ceed_data, X_loc, op_error, op_error_ctx)); 403 PetscScalar l2_error; 404 PetscCall(ComputeL2Error(X, &l2_error, op_error_ctx)); 405 406 if (!test_mode || l2_error > tolerance) { 407 PetscCall(MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm)); 408 PetscCall(MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm)); 409 PetscCall(PetscPrintf(comm, 410 " L2 Error : %e\n" 411 " CG Solve Time : %g (%g) sec\n", 412 (double)l2_error, rt_max, rt_min)); 413 } 414 } 415 if (benchmark_mode && (!test_mode)) { 416 PetscCall(PetscPrintf(comm, " DoFs/Sec in CG : %g (%g) million\n", 1e-6 * g_size * its / rt_max, 417 1e-6 * g_size * its / rt_min)); 418 } 419 } 420 421 // Output solution 422 if (write_solution) { 423 PetscViewer vtk_viewer_soln; 424 425 PetscCall(PetscViewerCreate(comm, &vtk_viewer_soln)); 426 PetscCall(PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK)); 427 PetscCall(PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu")); 428 PetscCall(VecView(X, vtk_viewer_soln)); 429 PetscCall(PetscViewerDestroy(&vtk_viewer_soln)); 430 } 431 432 // Cleanup 433 PetscCall(VecDestroy(&X)); 434 PetscCall(VecDestroy(&X_loc)); 435 PetscCall(VecDestroy(&op_apply_ctx->Y_loc)); 436 PetscCall(VecDestroy(&op_error_ctx->Y_loc)); 437 PetscCall(MatDestroy(&mat_O)); 438 PetscCall(PetscFree(op_apply_ctx)); 439 PetscCall(PetscFree(op_error_ctx)); 440 PetscCall(CeedDataDestroy(0, ceed_data)); 441 PetscCall(DMDestroy(&dm_mesh)); 442 PetscCall(DMDestroy(&dm_swarm)); 443 444 PetscCall(VecDestroy(&rhs)); 445 PetscCall(VecDestroy(&target)); 446 PetscCall(KSPDestroy(&ksp)); 447 CeedOperatorDestroy(&op_error); 448 CeedDestroy(&ceed); 449 return PetscFinalize(); 450 } 451