1 // Copyright (c) 2017-2026, 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
main(int argc,char ** argv)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 = VECSTANDARD;
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 // Set up libCEED
150 CeedInit(ceed_resource, &ceed);
151 CeedMemType mem_type_backend;
152 CeedGetPreferredMemType(ceed, &mem_type_backend);
153
154 // Set background mesh vec_type
155 switch (mem_type_backend) {
156 case CEED_MEM_HOST:
157 vec_type = VECSTANDARD;
158 break;
159 case CEED_MEM_DEVICE: {
160 const char *resolved;
161
162 CeedGetResource(ceed, &resolved);
163 if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
164 else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
165 else vec_type = VECSTANDARD;
166 }
167 }
168
169 // Setup DM
170 if (read_mesh) {
171 PetscCall(DMPlexCreateFromFile(comm, filename, NULL, PETSC_TRUE, &dm_mesh));
172 } else {
173 PetscCall(DMCreate(comm, &dm_mesh));
174 PetscCall(DMSetType(dm_mesh, DMPLEX));
175 PetscCall(DMSetFromOptions(dm_mesh));
176
177 // -- Check for tensor product mesh
178 {
179 PetscBool is_simplex;
180
181 PetscCall(DMPlexIsSimplex(dm_mesh, &is_simplex));
182 PetscCheck(!is_simplex, comm, PETSC_ERR_USER, "Only tensor-product background meshes supported");
183 }
184 }
185 PetscCall(DMSetVecType(dm_mesh, vec_type));
186 PetscCall(DMSetFromOptions(dm_mesh));
187
188 PetscCall(DMGetDimension(dm_mesh, &dim));
189 PetscCall(SetupDMByDegree(dm_mesh, degree, q_extra, num_comp_u, dim, bp_options[bp_choice].enforce_bc));
190
191 // View mesh
192 PetscCall(DMViewFromOptions(dm_mesh, NULL, "-dm_view"));
193
194 // Create particle swarm
195 PetscCall(DMCreate(comm, &dm_swarm));
196 PetscCall(DMSetType(dm_swarm, DMSWARM));
197 PetscCall(DMSetDimension(dm_swarm, dim));
198 PetscCall(DMSwarmSetType(dm_swarm, DMSWARM_PIC));
199 PetscCall(DMSwarmSetCellDM(dm_swarm, dm_mesh));
200
201 // -- Swarm field
202 PetscCall(DMSwarmRegisterPetscDatatypeField(dm_swarm, DMSwarmPICField_u, num_comp_u, PETSC_SCALAR));
203 PetscCall(DMSwarmFinalizeFieldRegister(dm_swarm));
204 {
205 PetscInt c_start, c_end, num_cells_local;
206 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &c_start, &c_end));
207 num_cells_local = c_end - c_start;
208 PetscCall(DMSwarmSetLocalSizes(dm_swarm, num_cells_local * num_points_per_cell, 0));
209 }
210 PetscCall(DMSetFromOptions(dm_swarm));
211
212 // -- Set swarm point locations
213 PetscCall(DMSwarmInitalizePointLocations(dm_swarm, point_swarm_type, num_points, num_points_per_cell));
214 PetscCall(DMSwarmVectorDefineField(dm_swarm, DMSwarmPICField_u));
215
216 // -- Final particle swarm
217 PetscCall(PetscObjectSetName((PetscObject)dm_swarm, "Particle Swarm"));
218 PetscCall(DMViewFromOptions(dm_swarm, NULL, "-dm_swarm_view"));
219
220 // Create vectors
221 PetscCall(DMCreateGlobalVector(dm_mesh, &X));
222 PetscCall(VecGetLocalSize(X, &l_size));
223 PetscCall(VecGetSize(X, &g_size));
224 PetscCall(DMCreateLocalVector(dm_mesh, &X_loc));
225 PetscCall(VecGetSize(X_loc, &xl_size));
226 PetscCall(VecDuplicate(X, &rhs));
227
228 // Operator
229 PetscCall(PetscMalloc1(1, &op_apply_ctx));
230 PetscCall(PetscMalloc1(1, &op_error_ctx));
231 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, op_apply_ctx, &mat_O));
232 PetscCall(MatSetDM(mat_O, dm_mesh));
233 PetscCall(MatShellSetOperation(mat_O, MATOP_MULT, (void (*)(void))MatMult_Ceed));
234 PetscCall(MatShellSetOperation(mat_O, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiag));
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_DIAGONAL));
306 } else {
307 PetscCall(PCSetType(pc, PCNONE));
308 }
309 PetscCall(KSPSetType(ksp, KSPCG));
310 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL));
311 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
312 }
313 PetscCall(KSPSetFromOptions(ksp));
314 PetscCall(KSPSetOperators(ksp, mat_O, mat_O));
315
316 // First run, if benchmarking
317 if (benchmark_mode) {
318 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1));
319 my_rt_start = MPI_Wtime();
320 PetscCall(KSPSolve(ksp, rhs, X));
321 my_rt = MPI_Wtime() - my_rt_start;
322 PetscCall(MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, comm));
323 // Set maxits based on first iteration timing
324 if (my_rt > 0.02) {
325 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 5));
326 } else {
327 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 20));
328 }
329 }
330
331 // Timed solve
332 PetscCall(VecZeroEntries(X));
333 PetscCall(PetscBarrier((PetscObject)ksp));
334
335 // -- Performance logging
336 PetscCall(PetscLogStageRegister("Solve Stage", &solve_stage));
337 PetscCall(PetscLogStagePush(solve_stage));
338
339 // -- Solve
340 my_rt_start = MPI_Wtime();
341 PetscCall(KSPSolve(ksp, rhs, X));
342 my_rt = MPI_Wtime() - my_rt_start;
343
344 // -- Performance logging
345 PetscCall(PetscLogStagePop());
346
347 // Output results
348 {
349 KSPType ksp_type;
350 KSPConvergedReason reason;
351 PetscReal rnorm;
352 PetscInt its;
353 PetscCall(KSPGetType(ksp, &ksp_type));
354 PetscCall(KSPGetConvergedReason(ksp, &reason));
355 PetscCall(KSPGetIterationNumber(ksp, &its));
356 PetscCall(KSPGetResidualNorm(ksp, &rnorm));
357 if (!test_mode || reason < 0 || rnorm > 1e-8) {
358 PetscCall(PetscPrintf(comm,
359 " KSP:\n"
360 " KSP Type : %s\n"
361 " KSP Convergence : %s\n"
362 " Total KSP Iterations : %" PetscInt_FMT "\n"
363 " Final rnorm : %e\n",
364 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm));
365 }
366 if (!test_mode) {
367 PetscCall(PetscPrintf(comm, " Performance:\n"));
368 }
369
370 // View true solution at particles
371 if (write_true_solution_swarm) {
372 Vec u_swarm, u_swarm_old;
373 PetscCall(DMSwarmSortGetAccess(dm_swarm));
374 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm));
375 PetscCall(VecDuplicate(u_swarm, &u_swarm_old));
376 PetscCall(VecCopy(u_swarm, u_swarm_old));
377 PetscCall(VecCopy(target, u_swarm));
378 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm));
379 PetscCall(DMSwarmSortRestoreAccess(dm_swarm));
380
381 PetscCall(DMSwarmViewXDMF(dm_swarm, "swarm.xmf"));
382
383 PetscCall(DMSwarmSortGetAccess(dm_swarm));
384 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm));
385 PetscCall(VecCopy(u_swarm_old, u_swarm));
386 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_u, &u_swarm));
387 PetscCall(DMSwarmSortRestoreAccess(dm_swarm));
388 PetscCall(VecDestroy(&u_swarm_old));
389 }
390
391 // View solution at mesh points
392 PetscCall(VecViewFromOptions(X, NULL, "-solution_view"));
393
394 // Compute L2 Error
395 {
396 // Set up error operator context
397 PetscCall(SetupErrorOperatorCtx(comm, dm_mesh, ceed, ceed_data, X_loc, op_error, op_error_ctx));
398 PetscScalar l2_error;
399 PetscCall(ComputeL2Error(X, &l2_error, op_error_ctx));
400
401 if (!test_mode || l2_error > tolerance) {
402 PetscCall(MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, comm));
403 PetscCall(MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, comm));
404 PetscCall(PetscPrintf(comm,
405 " L2 Error : %e\n"
406 " CG Solve Time : %g (%g) sec\n",
407 (double)l2_error, rt_max, rt_min));
408 }
409 }
410 if (benchmark_mode && (!test_mode)) {
411 PetscCall(PetscPrintf(comm, " DoFs/Sec in CG : %g (%g) million\n", 1e-6 * g_size * its / rt_max,
412 1e-6 * g_size * its / rt_min));
413 }
414 }
415
416 // Output solution
417 if (write_solution) {
418 PetscViewer vtk_viewer_soln;
419
420 PetscCall(PetscViewerCreate(comm, &vtk_viewer_soln));
421 PetscCall(PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK));
422 PetscCall(PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"));
423 PetscCall(VecView(X, vtk_viewer_soln));
424 PetscCall(PetscViewerDestroy(&vtk_viewer_soln));
425 }
426
427 // Cleanup
428 PetscCall(VecDestroy(&X));
429 PetscCall(VecDestroy(&X_loc));
430 PetscCall(VecDestroy(&op_apply_ctx->Y_loc));
431 PetscCall(VecDestroy(&op_error_ctx->Y_loc));
432 PetscCall(MatDestroy(&mat_O));
433 PetscCall(PetscFree(op_apply_ctx));
434 PetscCall(PetscFree(op_error_ctx));
435 PetscCall(CeedDataDestroy(0, ceed_data));
436 PetscCall(DMDestroy(&dm_mesh));
437 PetscCall(DMDestroy(&dm_swarm));
438
439 PetscCall(VecDestroy(&rhs));
440 PetscCall(VecDestroy(&target));
441 PetscCall(KSPDestroy(&ksp));
442 CeedOperatorDestroy(&op_error);
443 CeedDestroy(&ceed);
444 return PetscFinalize();
445 }
446