1 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
2 This example tests PCVPBJacobiSetBlocks().\n\n";
3
4 /*
5 Include "petscsnes.h" so that we can use SNES solvers. Note that this
6 file automatically includes:
7 petscsys.h - base PETSc routines petscvec.h - vectors
8 petscmat.h - matrices
9 petscis.h - index sets petscksp.h - Krylov subspace methods
10 petscviewer.h - viewers petscpc.h - preconditioners
11 petscksp.h - linear solvers
12 */
13
14 #include <petscsnes.h>
15
16 /*
17 User-defined routines
18 */
19 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
20 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
21 extern PetscErrorCode FormInitialGuess(Vec);
22
main(int argc,char ** argv)23 int main(int argc, char **argv)
24 {
25 SNES snes; /* SNES context */
26 Vec x, r, F, U; /* vectors */
27 Mat J; /* Jacobian matrix */
28 PetscInt its, n = 5, nb, maxit, maxf, *lens;
29 PetscMPIInt size;
30 PetscScalar h, xp, v, none = -1.0;
31 PetscReal abstol, rtol, stol, norm;
32 KSP ksp;
33 PC pc;
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
38 PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
39 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
40 PetscCheck(n % 5 == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "n must be a multiple of 5");
41 h = 1.0 / (n - 1);
42
43 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44 Create vector data structures
45 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46 /*
47 Note that we form 1 vector from scratch and then duplicate as needed.
48 */
49 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
50 PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
51 PetscCall(VecSetFromOptions(x));
52 PetscCall(VecDuplicate(x, &r));
53 PetscCall(VecDuplicate(x, &F));
54 PetscCall(VecDuplicate(x, &U));
55
56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 Create matrix data structure
58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59
60 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
61 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
62 PetscCall(MatSetFromOptions(J));
63 PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
64 PetscCall(MatSetBlockSize(J, 5));
65
66 nb = 3 * n / 5;
67 PetscCall(PetscMalloc1(nb, &lens));
68 for (PetscInt i = 0; i < nb / 3; i++) {
69 lens[3 * i + 0] = 1;
70 lens[3 * i + 1] = 2;
71 lens[3 * i + 2] = 2;
72 }
73 PetscCall(MatSetVariableBlockSizes(J, nb, lens));
74 PetscCall(PetscFree(lens));
75
76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 Create nonlinear solver context
78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79
80 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
81 PetscCall(SNESGetKSP(snes, &ksp));
82 PetscCall(KSPGetPC(ksp, &pc));
83 PetscCall(PCSetType(pc, PCVPBJACOBI));
84
85 /*
86 Set function evaluation routine and vector
87 */
88 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
89
90 /*
91 Set Jacobian matrix data structure and default Jacobian evaluation
92 routine. User can override with:
93 -snes_fd : default finite differencing approximation of Jacobian
94 -snes_mf : matrix-free Newton-Krylov method with no preconditioning
95 (unless user explicitly sets preconditioner)
96 -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
97 but use matrix-free approx for Jacobian-vector
98 products within Newton-Krylov method
99 */
100
101 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
102
103 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 Customize nonlinear solver; set runtime options
105 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106
107 /*
108 Set names for some vectors to facilitate monitoring (optional)
109 */
110 PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
111 PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
112
113 /*
114 Set SNES/KSP/KSP/PC runtime options, e.g.,
115 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
116 */
117 PetscCall(SNESSetFromOptions(snes));
118
119 /*
120 Print parameters used for convergence testing (optional) ... just
121 to demonstrate this routine; this information is also printed with
122 the option -snes_view
123 */
124 PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Initialize application:
129 Store right-hand side of PDE and exact solution
130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131
132 xp = 0.0;
133 for (PetscInt i = 0; i < n; i++) {
134 v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
135 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
136 v = xp * xp * xp;
137 PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
138 xp += h;
139 }
140
141 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 Evaluate initial guess; then solve nonlinear system
143 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 /*
145 Note: The user should initialize the vector, x, with the initial guess
146 for the nonlinear solver prior to calling SNESSolve(). In particular,
147 to employ an initial guess of zero, the user should explicitly set
148 this vector to zero by calling VecSet().
149 */
150 PetscCall(FormInitialGuess(x));
151 PetscCall(SNESSolve(snes, NULL, x));
152 PetscCall(SNESGetIterationNumber(snes, &its));
153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
154
155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156 Check solution and clean up
157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158
159 /*
160 Check the error
161 */
162 PetscCall(VecAXPY(x, none, U));
163 PetscCall(VecNorm(x, NORM_2, &norm));
164 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
165
166 /*
167 Free work space. All PETSc objects should be destroyed when they
168 are no longer needed.
169 */
170 PetscCall(VecDestroy(&x));
171 PetscCall(VecDestroy(&r));
172 PetscCall(VecDestroy(&U));
173 PetscCall(VecDestroy(&F));
174 PetscCall(MatDestroy(&J));
175 PetscCall(SNESDestroy(&snes));
176 PetscCall(PetscFinalize());
177 return 0;
178 }
179
180 /*
181 FormInitialGuess - Computes initial guess.
182
183 Input/Output Parameter:
184 . x - the solution vector
185 */
FormInitialGuess(Vec x)186 PetscErrorCode FormInitialGuess(Vec x)
187 {
188 PetscScalar pfive = .50;
189
190 PetscFunctionBeginUser;
191 PetscCall(VecSet(x, pfive));
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
195 /*
196 FormFunction - Evaluates nonlinear function, F(x).
197
198 Input Parameters:
199 . snes - the SNES context
200 . x - input vector
201 . ctx - optional user-defined context, as set by SNESSetFunction()
202
203 Output Parameter:
204 . f - function vector
205
206 Note:
207 The user-defined context can contain any application-specific data
208 needed for the function evaluation (such as various parameters, work
209 vectors, and grid information). In this program the context is just
210 a vector containing the right-hand side of the discretized PDE.
211 */
212
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)213 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
214 {
215 Vec g = (Vec)ctx;
216 const PetscScalar *xx, *gg;
217 PetscScalar *ff, d;
218 PetscInt i, n;
219
220 PetscFunctionBeginUser;
221 /*
222 Get pointers to vector data.
223 - For default PETSc vectors, VecGetArray() returns a pointer to
224 the data array. Otherwise, the routine is implementation dependent.
225 - You MUST call VecRestoreArray() when you no longer need access to
226 the array.
227 */
228 PetscCall(VecGetArrayRead(x, &xx));
229 PetscCall(VecGetArray(f, &ff));
230 PetscCall(VecGetArrayRead(g, &gg));
231
232 /*
233 Compute function
234 */
235 PetscCall(VecGetSize(x, &n));
236 d = (PetscReal)(n - 1);
237 d = d * d;
238 ff[0] = xx[0];
239 for (i = 1; i < n - 1; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - gg[i];
240 ff[n - 1] = xx[n - 1] - 1.0;
241
242 /*
243 Restore vectors
244 */
245 PetscCall(VecRestoreArrayRead(x, &xx));
246 PetscCall(VecRestoreArray(f, &ff));
247 PetscCall(VecRestoreArrayRead(g, &gg));
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
251 /*
252 FormJacobian - Evaluates Jacobian matrix.
253
254 Input Parameters:
255 . snes - the SNES context
256 . x - input vector
257 . dummy - optional user-defined context (not used here)
258
259 Output Parameters:
260 . jac - Jacobian matrix
261 . B - optionally different matrix used to construct the preconditioner
262
263 */
264
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)265 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
266 {
267 const PetscScalar *xx;
268 PetscScalar A[3], d;
269 PetscInt i, n, j[3];
270
271 PetscFunctionBeginUser;
272 /*
273 Get pointer to vector data
274 */
275 PetscCall(VecGetArrayRead(x, &xx));
276
277 /*
278 Compute Jacobian entries and insert into matrix.
279 - Note that in this case we set all elements for a particular
280 row at once.
281 */
282 PetscCall(VecGetSize(x, &n));
283 d = (PetscReal)(n - 1);
284 d = d * d;
285
286 /*
287 Interior grid points
288 */
289 for (i = 1; i < n - 1; i++) {
290 j[0] = i - 1;
291 j[1] = i;
292 j[2] = i + 1;
293 A[0] = d;
294 A[1] = -2.0 * d + 2.0 * xx[i];
295 A[2] = d;
296 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
297 }
298
299 /*
300 Boundary points
301 */
302 i = 0;
303 A[0] = 1.0;
304
305 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
306
307 i = n - 1;
308 A[0] = 1.0;
309
310 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
311
312 /*
313 Restore vector
314 */
315 PetscCall(VecRestoreArrayRead(x, &xx));
316
317 /*
318 Assemble matrix
319 */
320 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
321 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
322 if (jac != B) {
323 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
324 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
325 }
326 PetscFunctionReturn(PETSC_SUCCESS);
327 }
328
329 /*TEST
330
331 testset:
332 args: -snes_monitor_short -snes_view -ksp_monitor
333 output_file: output/ex5_1.out
334 filter: grep -v "type: seqaij"
335
336 test:
337 suffix: 1
338
339 test:
340 suffix: cuda
341 requires: cuda
342 args: -mat_type aijcusparse -vec_type cuda
343
344 test:
345 suffix: kok
346 requires: kokkos_kernels
347 args: -mat_type aijkokkos -vec_type kokkos
348
349 # this is just a test for SNESKSPTRANSPOSEONLY and KSPSolveTranspose() to behave properly
350 # the solution is wrong on purpose
351 test:
352 requires: !single !complex
353 suffix: transpose_only
354 args: -snes_monitor_short -snes_view -ksp_monitor -snes_type ksptransposeonly -pc_type ilu -snes_test_jacobian -snes_test_jacobian_view -ksp_view_rhs -ksp_view_solution -ksp_view_mat_explicit -ksp_view_preconditioned_operator_explicit
355
356 test:
357 requires: mumps
358 suffix: mumps
359 args: -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_15 1 -snes_monitor_short -ksp_monitor
360
361 test:
362 suffix: fieldsplit_1
363 args: -snes_monitor_short -snes_view -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 2,3,4
364
365 test:
366 suffix: fieldsplit_2
367 args: -snes_monitor_short -snes_view -ksp_monitor -pc_type fieldsplit -pc_fieldsplit_0_fields 1,0 -pc_fieldsplit_1_fields 2,3,4
368
369 TEST*/
370