Actual source code: ex13.c
1: static char help[] = "Benchmark Poisson Problem in 2d and 3d with finite elements.\n\
2: We solve the Poisson problem in a rectangular domain\n\
3: using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: #include <petscdmplex.h>
6: #include <petscsnes.h>
7: #include <petscds.h>
8: #include <petscconvest.h>
9: #if PetscDefined(HAVE_AMGX)
10: #include <amgx_c.h>
11: #endif
13: typedef struct {
14: PetscInt nit; /* Number of benchmark iterations */
15: PetscBool strong; /* Do not integrate the Laplacian by parts */
16: } AppCtx;
18: static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
19: {
20: PetscInt d;
21: *u = 0.0;
22: for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
23: return PETSC_SUCCESS;
24: }
26: static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
27: {
28: PetscInt d;
29: for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
30: }
32: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
33: {
34: for (PetscInt d = 0; d < dim; ++d) f1[d] = u_x[d];
35: }
37: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
38: {
39: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
40: }
42: static PetscErrorCode quadratic_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
43: {
44: *u = PetscSqr(x[0]) + PetscSqr(x[1]);
45: return PETSC_SUCCESS;
46: }
48: static void f0_strong_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
49: {
50: for (PetscInt d = 0; d < dim; ++d) f0[0] -= u_x[dim + d * dim + d];
51: f0[0] += 4.0;
52: }
54: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
55: {
56: PetscFunctionBeginUser;
57: options->nit = 10;
58: options->strong = PETSC_FALSE;
59: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
60: PetscCall(PetscOptionsInt("-benchmark_it", "Solve the benchmark problem this many times", "ex13.c", options->nit, &options->nit, NULL));
61: PetscCall(PetscOptionsBool("-strong", "Do not integrate the Laplacian by parts", "ex13.c", options->strong, &options->strong, NULL));
62: PetscOptionsEnd();
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
67: {
68: PetscFunctionBeginUser;
69: PetscCall(DMCreate(comm, dm));
70: PetscCall(DMSetType(*dm, DMPLEX));
71: PetscCall(DMSetFromOptions(*dm));
72: PetscCall(DMSetApplicationContext(*dm, user));
73: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
74: { // perturb to get general coordinates
75: Vec coordinates;
76: PetscScalar *coords;
77: PetscInt nloc;
78: PetscRandom rnd;
79: PetscReal del;
80: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd));
81: PetscCall(PetscRandomSetInterval(rnd, -PETSC_SQRT_MACHINE_EPSILON, PETSC_SQRT_MACHINE_EPSILON));
82: PetscCall(PetscRandomSetFromOptions(rnd));
83: PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
84: PetscCall(VecGetArray(coordinates, &coords));
85: PetscCall(VecGetLocalSize(coordinates, &nloc));
86: for (PetscInt v = 0; v < nloc; ++v) {
87: PetscCall(PetscRandomGetValueReal(rnd, &del));
88: coords[v] += del * coords[v];
89: }
90: PetscCall(VecRestoreArray(coordinates, &coords));
91: PetscCall(PetscRandomDestroy(&rnd));
92: }
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
97: {
98: PetscDS ds;
99: DMLabel label;
100: const PetscInt id = 1;
102: PetscFunctionBeginUser;
103: PetscCall(DMGetDS(dm, &ds));
104: PetscCall(DMGetLabel(dm, "marker", &label));
105: if (user->strong) {
106: PetscCall(PetscDSSetResidual(ds, 0, f0_strong_u, NULL));
107: PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u, user));
108: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscFortranCallbackFn *)quadratic_u, NULL, user, NULL));
109: } else {
110: PetscCall(PetscDSSetResidual(ds, 0, f0_trig_u, f1_u));
111: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu));
112: PetscCall(PetscDSSetExactSolution(ds, 0, trig_u, user));
113: PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)trig_u, NULL, user, NULL));
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
119: {
120: DM cdm = dm;
121: PetscFE fe;
122: DMPolytopeType ct;
123: PetscBool simplex;
124: PetscInt dim, cStart;
125: char prefix[PETSC_MAX_PATH_LEN];
127: PetscFunctionBeginUser;
128: PetscCall(DMGetDimension(dm, &dim));
129: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
130: PetscCall(DMPlexGetCellType(dm, cStart, &ct));
131: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; // false
132: /* Create finite element */
133: PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name));
134: PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe));
135: PetscCall(PetscObjectSetName((PetscObject)fe, name));
136: /* Set discretization and boundary conditions for each mesh */
137: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
138: PetscCall(DMCreateDS(dm));
139: PetscCall((*setup)(dm, user));
140: while (cdm) {
141: PetscCall(DMCopyDisc(dm, cdm));
142: /* TODO: Check whether the boundary of coarse meshes is marked */
143: PetscCall(DMGetCoarseDM(cdm, &cdm));
144: }
145: PetscCall(PetscFEDestroy(&fe));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: int main(int argc, char **argv)
150: {
151: DM dm; /* Problem specification */
152: SNES snes; /* Nonlinear solver */
153: Vec u; /* Solutions */
154: AppCtx user; /* User-defined work context */
155: PetscLogDouble time;
156: Mat Amat;
158: PetscFunctionBeginUser;
159: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
160: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
161: /* system */
162: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
163: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
164: PetscCall(SNESSetDM(snes, dm));
165: PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user));
166: PetscCall(DMCreateGlobalVector(dm, &u));
167: {
168: PetscInt N;
169: PetscCall(VecGetSize(u, &N));
170: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number equations N = %" PetscInt_FMT "\n", N));
171: }
172: PetscCall(SNESSetFromOptions(snes));
173: PetscCall(PetscObjectSetName((PetscObject)u, "potential"));
174: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
175: PetscCall(DMSNESCheckFromOptions(snes, u));
176: PetscCall(PetscTime(&time));
177: PetscCall(SNESSetUp(snes));
178: #if PetscDefined(HAVE_AMGX)
179: KSP ksp;
180: PC pc;
181: PetscBool flg;
182: AMGX_resources_handle rsc;
183: PetscCall(SNESGetKSP(snes, &ksp));
184: PetscCall(KSPGetPC(ksp, &pc));
185: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCAMGX, &flg));
186: if (flg) {
187: PetscCall(PCAmgXGetResources(pc, (void *)&rsc));
188: /* do ... with resource */
189: }
190: #endif
191: PetscCall(SNESGetJacobian(snes, &Amat, NULL, NULL, NULL));
192: PetscCall(MatSetOption(Amat, MAT_SPD, PETSC_TRUE));
193: PetscCall(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE));
194: PetscCall(SNESSolve(snes, NULL, u));
195: PetscCall(PetscTimeSubtract(&time));
196: /* Benchmark system */
197: if (user.nit) {
198: Vec b;
199: PetscLogStage kspstage;
200: PetscCall(PetscLogStageRegister("Solve only", &kspstage));
201: PetscCall(PetscLogStagePush(kspstage));
202: PetscCall(SNESGetSolution(snes, &u));
203: PetscCall(SNESGetFunction(snes, &b, NULL, NULL));
204: for (PetscInt i = 0; i < user.nit; i++) {
205: PetscCall(VecZeroEntries(u));
206: PetscCall(SNESSolve(snes, NULL, u));
207: }
208: PetscCall(PetscLogStagePop());
209: }
210: PetscCall(SNESGetSolution(snes, &u));
211: PetscCall(VecViewFromOptions(u, NULL, "-potential_view"));
212: /* Cleanup */
213: PetscCall(VecDestroy(&u));
214: PetscCall(SNESDestroy(&snes));
215: PetscCall(DMDestroy(&dm));
216: PetscCall(PetscFinalize());
217: return 0;
218: }
220: /*TEST
222: test:
223: suffix: strong
224: requires: triangle
225: args: -dm_plex_dim 2 -dm_refine 1 -benchmark_it 0 -dmsnes_check -potential_petscspace_degree 2 -dm_ds_jet_degree 2 -strong -pc_type jacobi
227: testset:
228: nsize: 4
229: output_file: output/ex13_comparison.out
230: args: -dm_plex_dim 3 -benchmark_it 2 -dm_plex_simplex 0 -dm_plex_box_faces 2,2,1 -dm_refine 2 -petscpartitioner_simple_node_grid 1,1,1 -petscpartitioner_simple_process_grid 2,2,1 -potential_petscspace_degree 2 -petscpartitioner_type simple -snes_type ksponly -dm_view -ksp_type cg -ksp_rtol 1e-12 -snes_lag_jacobian -2 -dm_plex_box_upper 2,2,1 -dm_plex_box_lower 0,0,0 -pc_type gamg -pc_gamg_process_eq_limit 200 -pc_gamg_coarse_eq_limit 1000 -pc_gamg_esteig_ksp_type cg -mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.05 -pc_gamg_reuse_interpolation true -pc_gamg_aggressive_square_graph true -pc_gamg_threshold 0.04 -pc_gamg_threshold_scale .25 -pc_gamg_aggressive_coarsening 2 -pc_gamg_mis_k_minimum_degree_ordering true -ksp_monitor -ksp_norm_type unpreconditioned
231: test:
232: suffix: comparison
233: test:
234: suffix: cuda
235: requires: cuda
236: args: -dm_mat_type aijcusparse -dm_vec_type cuda
237: test:
238: suffix: kokkos
239: requires: kokkos_kernels
240: args: -dm_mat_type aijkokkos -dm_vec_type kokkos
241: test:
242: suffix: kokkos_sycl
243: requires: sycl kokkos_kernels
244: args: -dm_mat_type aijkokkos -dm_vec_type kokkos
245: test:
246: suffix: aijmkl_comp
247: requires: mkl_sparse
248: args: -dm_mat_type aijmkl
250: testset:
251: requires: cuda amgx
252: filter: grep -v Built | grep -v "AMGX version" | grep -v "CUDA Runtime"
253: output_file: output/ex13_amgx.out
254: args: -dm_plex_dim 2 -dm_plex_box_faces 2,2 -dm_refine 2 -petscpartitioner_type simple -potential_petscspace_degree 2 -dm_plex_simplex 0 -ksp_monitor \
255: -snes_type ksponly -dm_view -ksp_type cg -ksp_norm_type unpreconditioned -ksp_converged_reason -snes_rtol 1.e-4 -pc_type amgx -benchmark_it 1 -pc_amgx_verbose false
256: nsize: 4
257: test:
258: suffix: amgx
259: args: -dm_mat_type aijcusparse -dm_vec_type cuda
260: test:
261: suffix: amgx_cpu
262: args: -dm_mat_type aij
264: TEST*/