xref: /petsc/src/snes/tests/ex1.c (revision 4e8208cbcbc709572b8abe32f33c78b69c819375) !
1 static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\
2 This example also illustrates the use of matrix coloring.  Runtime options include:\n\
3   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
4      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
5   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
6   -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
7 
8 /*
9 
10     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
11     the partial differential equation
12 
13             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
14 
15     with boundary conditions
16 
17              u = 0  for  x = 0, x = 1, y = 0, y = 1.
18 
19     A finite difference approximation with the usual 5-point stencil
20     is used to discretize the boundary value problem to obtain a nonlinear
21     system of equations.
22 
23     The parallel version of this code is snes/tutorials/ex5.c
24 
25 */
26 
27 /*
28    Include "petscsnes.h" so that we can use SNES solvers.  Note that
29    this file automatically includes:
30      petscsys.h       - base PETSc routines   petscvec.h - vectors
31      petscmat.h - matrices
32      petscis.h     - index sets            petscksp.h - Krylov subspace methods
33      petscviewer.h - viewers               petscpc.h  - preconditioners
34      petscksp.h   - linear solvers
35 */
36 
37 #include <petscsnes.h>
38 
39 /*
40    User-defined application context - contains data needed by the
41    application-provided call-back routines, FormJacobian() and
42    FormFunction().
43 */
44 typedef struct {
45   PetscReal param; /* test problem parameter */
46   PetscInt  mx;    /* Discretization in x-direction */
47   PetscInt  my;    /* Discretization in y-direction */
48 } AppCtx;
49 
50 /*
51    User-defined routines
52 */
53 static PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
54 static PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
55 static PetscErrorCode FormInitialGuess(AppCtx *, Vec);
56 static PetscErrorCode ConvergenceTest(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
57 static PetscErrorCode ConvergenceDestroy(PetscCtxRt);
58 static PetscErrorCode postcheck(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
59 static PetscErrorCode monitor_change_deltamax(SNES, PetscInt, PetscReal, void *);
60 
main(int argc,char ** argv)61 int main(int argc, char **argv)
62 {
63   SNES          snes; /* nonlinear solver context */
64   Vec           x, r; /* solution, residual vectors */
65   Mat           J;    /* Jacobian matrix */
66   AppCtx        user; /* user-defined application context */
67   PetscInt      i, its, N, hist_its[50];
68   PetscMPIInt   size;
69   PetscReal     bratu_lambda_max = 6.81, bratu_lambda_min = 0., history[50];
70   MatFDColoring fdcoloring;
71   PetscBool     matrix_free = PETSC_FALSE, flg, fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE, pc = PETSC_FALSE, prunejacobian = PETSC_FALSE, null_appctx = PETSC_TRUE, test_tr_deltamax = PETSC_FALSE;
72   KSP           ksp;
73   PetscInt     *testarray;
74 
75   PetscFunctionBeginUser;
76   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
77   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
78   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
79 
80   /*
81      Initialize problem parameters
82   */
83   user.mx    = 4;
84   user.my    = 4;
85   user.param = 6.0;
86   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
87   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
88   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL));
89   PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL));
90   PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range");
91   N = user.mx * user.my;
92   PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL));
93   PetscCall(PetscOptionsGetBool(NULL, NULL, "-prune_jacobian", &prunejacobian, NULL));
94   PetscCall(PetscOptionsGetBool(NULL, NULL, "-null_appctx", &null_appctx, NULL));
95   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_tr_deltamax", &test_tr_deltamax, NULL));
96 
97   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98      Create nonlinear solver context
99      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 
101   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
102   if (pc) {
103     PetscCall(SNESSetType(snes, SNESNEWTONTR));
104     PetscCall(SNESNewtonTRSetPostCheck(snes, postcheck, NULL));
105   }
106   if (test_tr_deltamax) {
107     PetscCall(SNESSetType(snes, SNESNEWTONTR));
108     PetscCall(SNESMonitorSet(snes, monitor_change_deltamax, NULL, NULL));
109   }
110 
111   /* Test application context handling from Python */
112   if (!null_appctx) PetscCall(SNESSetApplicationContext(snes, (void *)&user));
113 
114   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115      Create vector data structures; set function evaluation routine
116      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117 
118   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
119   PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
120   PetscCall(VecSetFromOptions(x));
121   PetscCall(VecDuplicate(x, &r));
122 
123   /*
124      Set function evaluation routine and vector.  Whenever the nonlinear
125      solver needs to evaluate the nonlinear function, it will call this
126      routine.
127       - Note that the final routine argument is the user-defined
128         context that provides application-specific data for the
129         function evaluation routine.
130   */
131   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));
132 
133   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134      Create matrix data structure; set Jacobian evaluation routine
135      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 
137   /*
138      Create matrix. Here we only approximately preallocate storage space
139      for the Jacobian.  See the users manual for a discussion of better
140      techniques for preallocating matrix memory.
141   */
142   PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL));
143   if (!matrix_free) {
144     PetscBool matrix_free_operator = PETSC_FALSE;
145     PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL));
146     if (matrix_free_operator) matrix_free = PETSC_FALSE;
147   }
148   if (!matrix_free) PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J));
149 
150   /*
151      This option will cause the Jacobian to be computed via finite differences
152     efficiently using a coloring of the columns of the matrix.
153   */
154   PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL));
155   PetscCheck(!matrix_free || !fd_coloring, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Use only one of -snes_mf, -snes_fd_coloring options! You can do -snes_mf_operator -snes_fd_coloring");
156 
157   if (fd_coloring) {
158     ISColoring  iscoloring;
159     MatColoring mc;
160     if (prunejacobian) {
161       /* Initialize x with random nonzero values so that the nonzeros in the Jacobian
162          can better reflect the sparsity structure of the Jacobian. */
163       PetscRandom rctx;
164       PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
165       PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0));
166       PetscCall(VecSetRandom(x, rctx));
167       PetscCall(PetscRandomDestroy(&rctx));
168     }
169 
170     /*
171       This initializes the nonzero structure of the Jacobian. This is artificial
172       because clearly if we had a routine to compute the Jacobian we won't need
173       to use finite differences.
174     */
175     PetscCall(FormJacobian(snes, x, J, J, &user));
176 
177     /*
178        Color the matrix, i.e. determine groups of columns that share no common
179       rows. These columns in the Jacobian can all be computed simultaneously.
180     */
181     PetscCall(MatColoringCreate(J, &mc));
182     PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
183     PetscCall(MatColoringSetFromOptions(mc));
184     PetscCall(MatColoringApply(mc, &iscoloring));
185     PetscCall(MatColoringDestroy(&mc));
186     /*
187        Create the data structure that SNESComputeJacobianDefaultColor() uses
188        to compute the actual Jacobians via finite differences.
189     */
190     PetscCall(MatFDColoringCreate(J, iscoloring, &fdcoloring));
191     PetscCall(MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn *)FormFunction, &user));
192     PetscCall(MatFDColoringSetFromOptions(fdcoloring));
193     PetscCall(MatFDColoringSetUp(J, iscoloring, fdcoloring));
194     /*
195         Tell SNES to use the routine SNESComputeJacobianDefaultColor()
196       to compute Jacobians.
197     */
198     PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring));
199     PetscCall(ISColoringDestroy(&iscoloring));
200     if (prunejacobian) PetscCall(SNESPruneJacobianColor(snes, J, J));
201   }
202   /*
203      Set Jacobian matrix data structure and default Jacobian evaluation
204      routine.  Whenever the nonlinear solver needs to compute the
205      Jacobian matrix, it will call this routine.
206       - Note that the final routine argument is the user-defined
207         context that provides application-specific data for the
208         Jacobian evaluation routine.
209       - The user can override with:
210          -snes_fd : default finite differencing approximation of Jacobian
211          -snes_mf : matrix-free Newton-Krylov method with no preconditioning
212                     (unless user explicitly sets preconditioner)
213          -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
214                              but use matrix-free approx for Jacobian-vector
215                              products within Newton-Krylov method
216   */
217   else if (!matrix_free) {
218     PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user));
219   }
220 
221   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222      Customize nonlinear solver; set runtime options
223    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224 
225   /*
226      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
227   */
228   PetscCall(SNESSetFromOptions(snes));
229 
230   /*
231      Set array that saves the function norms.  This array is intended
232      when the user wants to save the convergence history for later use
233      rather than just to view the function norms via -snes_monitor.
234   */
235   PetscCall(SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE));
236 
237   /*
238       Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
239       user provided test before the specialized test. The convergence context is just an array to
240       test that it gets properly freed at the end
241   */
242   if (use_convergence_test) {
243     PetscCall(SNESGetKSP(snes, &ksp));
244     PetscCall(PetscMalloc1(5, &testarray));
245     PetscCall(KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy));
246   }
247 
248   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249      Evaluate initial guess; then solve nonlinear system
250    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
251   /*
252      Note: The user should initialize the vector, x, with the initial guess
253      for the nonlinear solver prior to calling SNESSolve().  In particular,
254      to employ an initial guess of zero, the user should explicitly set
255      this vector to zero by calling VecSet().
256   */
257   PetscCall(FormInitialGuess(&user, x));
258   PetscCall(SNESSolve(snes, NULL, x));
259   PetscCall(SNESGetIterationNumber(snes, &its));
260   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
261 
262   /*
263      Print the convergence history.  This is intended just to demonstrate
264      use of the data attained via SNESSetConvergenceHistory().
265   */
266   PetscCall(PetscOptionsHasName(NULL, NULL, "-print_history", &flg));
267   if (flg) {
268     for (i = 0; i < its + 1; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iteration %" PetscInt_FMT ": Linear iterations %" PetscInt_FMT " Function norm = %g\n", i, hist_its[i], (double)history[i]));
269   }
270 
271   /* Test NewtonTR API */
272   PetscCall(SNESNewtonTRSetTolerances(snes, 1.0, 2.0, 3.0));
273   PetscCall(SNESNewtonTRSetUpdateParameters(snes, 4.0, 5.0, 6.0, 7.0, 8.0));
274   PetscCall(PetscObjectTypeCompare((PetscObject)snes, SNESNEWTONTR, &flg));
275   if (flg) {
276     PetscReal tmp[5];
277 
278     PetscCall(SNESNewtonTRGetTolerances(snes, &tmp[0], &tmp[1], &tmp[2]));
279     PetscCheck(tmp[0] == 1.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
280     PetscCheck(tmp[1] == 2.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
281     PetscCheck(tmp[2] == 3.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
282     PetscCall(SNESNewtonTRGetUpdateParameters(snes, &tmp[0], &tmp[1], &tmp[2], &tmp[3], &tmp[4]));
283     PetscCheck(tmp[0] == 4.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
284     PetscCheck(tmp[1] == 5.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
285     PetscCheck(tmp[2] == 6.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
286     PetscCheck(tmp[3] == 7.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
287     PetscCheck(tmp[4] == 8.0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Wrong value");
288   }
289 
290   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291      Free work space.  All PETSc objects should be destroyed when they
292      are no longer needed.
293    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
294 
295   if (!matrix_free) PetscCall(MatDestroy(&J));
296   if (fd_coloring) PetscCall(MatFDColoringDestroy(&fdcoloring));
297   PetscCall(VecDestroy(&x));
298   PetscCall(VecDestroy(&r));
299   PetscCall(SNESDestroy(&snes));
300   PetscCall(PetscFinalize());
301   return 0;
302 }
303 
304 /*
305    FormInitialGuess - Forms initial approximation.
306 
307    Input Parameters:
308    user - user-defined application context
309    X - vector
310 
311    Output Parameter:
312    X - vector
313  */
FormInitialGuess(AppCtx * user,Vec X)314 PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
315 {
316   PetscInt     i, j, row, mx, my;
317   PetscReal    lambda, temp1, temp, hx, hy;
318   PetscScalar *x;
319 
320   PetscFunctionBeginUser;
321   mx     = user->mx;
322   my     = user->my;
323   lambda = user->param;
324 
325   hx = 1.0 / (PetscReal)(mx - 1);
326   hy = 1.0 / (PetscReal)(my - 1);
327 
328   /*
329      Get a pointer to vector data.
330        - For default PETSc vectors, VecGetArray() returns a pointer to
331          the data array.  Otherwise, the routine is implementation dependent.
332        - You MUST call VecRestoreArray() when you no longer need access to
333          the array.
334   */
335   PetscCall(VecGetArray(X, &x));
336   temp1 = lambda / (lambda + 1.0);
337   for (j = 0; j < my; j++) {
338     temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
339     for (i = 0; i < mx; i++) {
340       row = i + j * mx;
341       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
342         x[row] = 0.0;
343         continue;
344       }
345       x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
346     }
347   }
348 
349   /*
350      Restore vector
351   */
352   PetscCall(VecRestoreArray(X, &x));
353   PetscFunctionReturn(PETSC_SUCCESS);
354 }
355 
356 /*
357    FormFunction - Evaluates nonlinear function, F(x).
358 
359    Input Parameters:
360 .  snes - the SNES context
361 .  X - input vector
362 .  ptr - optional user-defined context, as set by SNESSetFunction()
363 
364    Output Parameter:
365 .  F - function vector
366  */
FormFunction(SNES snes,Vec X,Vec F,void * ptr)367 PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
368 {
369   AppCtx            *user = (AppCtx *)ptr;
370   PetscInt           i, j, row, mx, my;
371   PetscReal          two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx;
372   PetscScalar        ut, ub, ul, ur, u, uxx, uyy, sc, *f;
373   const PetscScalar *x;
374 
375   PetscFunctionBeginUser;
376   mx     = user->mx;
377   my     = user->my;
378   lambda = user->param;
379   hx     = one / (PetscReal)(mx - 1);
380   hy     = one / (PetscReal)(my - 1);
381   sc     = hx * hy;
382   hxdhy  = hx / hy;
383   hydhx  = hy / hx;
384 
385   /*
386      Get pointers to vector data
387   */
388   PetscCall(VecGetArrayRead(X, &x));
389   PetscCall(VecGetArray(F, &f));
390 
391   /*
392      Compute function
393   */
394   for (j = 0; j < my; j++) {
395     for (i = 0; i < mx; i++) {
396       row = i + j * mx;
397       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
398         f[row] = x[row];
399         continue;
400       }
401       u      = x[row];
402       ub     = x[row - mx];
403       ul     = x[row - 1];
404       ut     = x[row + mx];
405       ur     = x[row + 1];
406       uxx    = (-ur + two * u - ul) * hydhx;
407       uyy    = (-ut + two * u - ub) * hxdhy;
408       f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u);
409     }
410   }
411 
412   /*
413      Restore vectors
414   */
415   PetscCall(VecRestoreArrayRead(X, &x));
416   PetscCall(VecRestoreArray(F, &f));
417   PetscFunctionReturn(PETSC_SUCCESS);
418 }
419 
420 /*
421    FormJacobian - Evaluates Jacobian matrix.
422 
423    Input Parameters:
424 .  snes - the SNES context
425 .  x - input vector
426 .  ptr - optional user-defined context, as set by SNESSetJacobian()
427 
428    Output Parameters:
429 .  A - Jacobian matrix
430 .  B - optionally different matrix used to construct the preconditioner
431 
432 */
FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void * ptr)433 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
434 {
435   AppCtx            *user = (AppCtx *)ptr; /* user-defined application context */
436   PetscInt           i, j, row, mx, my, col[5];
437   PetscScalar        two = 2.0, one = 1.0, lambda, v[5], sc;
438   const PetscScalar *x;
439   PetscReal          hx, hy, hxdhy, hydhx;
440 
441   PetscFunctionBeginUser;
442   mx     = user->mx;
443   my     = user->my;
444   lambda = user->param;
445   hx     = 1.0 / (PetscReal)(mx - 1);
446   hy     = 1.0 / (PetscReal)(my - 1);
447   sc     = hx * hy;
448   hxdhy  = hx / hy;
449   hydhx  = hy / hx;
450 
451   /*
452      Get pointer to vector data
453   */
454   PetscCall(VecGetArrayRead(X, &x));
455 
456   /*
457      Compute entries of the Jacobian
458   */
459   for (j = 0; j < my; j++) {
460     for (i = 0; i < mx; i++) {
461       row = i + j * mx;
462       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
463         PetscCall(MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES));
464         continue;
465       }
466       v[0]   = -hxdhy;
467       col[0] = row - mx;
468       v[1]   = -hydhx;
469       col[1] = row - 1;
470       v[2]   = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]);
471       col[2] = row;
472       v[3]   = -hydhx;
473       col[3] = row + 1;
474       v[4]   = -hxdhy;
475       col[4] = row + mx;
476       PetscCall(MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES));
477     }
478   }
479 
480   /*
481      Restore vector
482   */
483   PetscCall(VecRestoreArrayRead(X, &x));
484 
485   /*
486      Assemble matrix
487   */
488   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
489   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
490 
491   if (jac != J) {
492     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
493     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
494   }
495   PetscFunctionReturn(PETSC_SUCCESS);
496 }
497 
ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason * reason,PetscCtx ctx)498 PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, PetscCtx ctx)
499 {
500   PetscFunctionBeginUser;
501   *reason = KSP_CONVERGED_ITERATING;
502   if (it > 1) {
503     *reason = KSP_CONVERGED_ITS;
504     PetscCall(PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n"));
505   }
506   PetscFunctionReturn(PETSC_SUCCESS);
507 }
508 
ConvergenceDestroy(PetscCtxRt ctx)509 PetscErrorCode ConvergenceDestroy(PetscCtxRt ctx)
510 {
511   PetscFunctionBeginUser;
512   PetscCall(PetscInfo(NULL, "User provided convergence destroy called\n"));
513   PetscCall(PetscFree(*(void **)ctx));
514   PetscFunctionReturn(PETSC_SUCCESS);
515 }
516 
postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool * changed_y,PetscBool * changed_w,PetscCtx ctx)517 PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, PetscCtx ctx)
518 {
519   PetscReal norm;
520   Vec       tmp;
521 
522   PetscFunctionBeginUser;
523   PetscCall(VecDuplicate(x, &tmp));
524   PetscCall(VecWAXPY(tmp, -1.0, x, w));
525   PetscCall(VecNorm(tmp, NORM_2, &norm));
526   PetscCall(VecDestroy(&tmp));
527   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm));
528   PetscFunctionReturn(PETSC_SUCCESS);
529 }
530 
monitor_change_deltamax(SNES snes,PetscInt it,PetscReal fnorm,PetscCtx ctx)531 PetscErrorCode monitor_change_deltamax(SNES snes, PetscInt it, PetscReal fnorm, PetscCtx ctx)
532 {
533   PetscFunctionBeginUser;
534   if (it == 0) PetscCall(SNESNewtonTRSetTolerances(snes, PETSC_CURRENT, 0.01, PETSC_CURRENT));
535   PetscFunctionReturn(PETSC_SUCCESS);
536 }
537 
538 /*TEST
539 
540    build:
541       requires: !single
542 
543    test:
544       args: -ksp_gmres_cgs_refinement_type refine_always
545 
546    test:
547       suffix: 2
548       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
549 
550    test:
551       suffix: 2_trdeltamax_change
552       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -pc -test_tr_deltamax
553 
554    test:
555       suffix: 2a
556       filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
557       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
558       requires: defined(PETSC_USE_INFO)
559 
560    test:
561       suffix: 2b
562       filter: grep -i  "User provided convergence test" > /dev/null  && echo "Found User provided convergence test"
563       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
564       requires: defined(PETSC_USE_INFO)
565 
566    test:
567       suffix: 2c
568       args: -snes_converged_reason -snes_type newtontr -snes_tr_qn {{same different}separate output} -pc_type mat -snes_view -snes_tr_qn_mat_type lmvmdfp -snes_tr_norm_type infinity
569 
570    test:
571       suffix: 3
572       args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
573 
574    test:
575       suffix: 4
576       args: -pc -par 6.807 -snes_monitor -snes_converged_reason
577 
578    test:
579       suffix: 5
580       args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always -prune_jacobian
581       output_file: output/ex1_3.out
582 
583    test:
584       suffix: 6
585       args: -snes_monitor draw:image:testfile -viewer_view
586 
587    test:
588       suffix: python
589       requires: petsc4py
590       args: -python -snes_type python -snes_python_type ex1.py:MySNES -snes_view -null_appctx {{0 1}separate output}
591       localrunfiles: ex1.py
592 
593 TEST*/
594