1 
2 #include <petsc/private/snesimpl.h>       /*I   "petsc/private/snesimpl.h"   I*/
3 #include <petscdm.h>
4 #include <petscsection.h>
5 #include <petscblaslapack.h>
6 
7 /*@C
8    SNESMonitorSolution - Monitors progress of the SNES solvers by calling
9    VecView() for the approximate solution at each iteration.
10 
11    Collective on SNES
12 
13    Input Parameters:
14 +  snes - the SNES context
15 .  its - iteration number
16 .  fgnorm - 2-norm of residual
17 -  dummy -  a viewer
18 
19    Options Database Keys:
20 .  -snes_monitor_solution [ascii binary draw][:filename][:viewer format] - plots solution at each iteration
21 
22    Level: intermediate
23 
24 .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
25 @*/
SNESMonitorSolution(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)26 PetscErrorCode  SNESMonitorSolution(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
27 {
28   PetscErrorCode ierr;
29   Vec            x;
30   PetscViewer    viewer = vf->viewer;
31 
32   PetscFunctionBegin;
33   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
34   ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr);
35   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
36   ierr = VecView(x,viewer);CHKERRQ(ierr);
37   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
38   PetscFunctionReturn(0);
39 }
40 
41 /*@C
42    SNESMonitorResidual - Monitors progress of the SNES solvers by calling
43    VecView() for the residual at each iteration.
44 
45    Collective on SNES
46 
47    Input Parameters:
48 +  snes - the SNES context
49 .  its - iteration number
50 .  fgnorm - 2-norm of residual
51 -  dummy -  a viewer
52 
53    Options Database Keys:
54 .  -snes_monitor_residual [ascii binary draw][:filename][:viewer format] - plots residual (not its norm) at each iteration
55 
56    Level: intermediate
57 
58 .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
59 @*/
SNESMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)60 PetscErrorCode  SNESMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
61 {
62   PetscErrorCode ierr;
63   Vec            x;
64   PetscViewer    viewer = vf->viewer;
65 
66   PetscFunctionBegin;
67   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
68   ierr = SNESGetFunction(snes,&x,NULL,NULL);CHKERRQ(ierr);
69   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
70   ierr = VecView(x,viewer);CHKERRQ(ierr);
71   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
72   PetscFunctionReturn(0);
73 }
74 
75 /*@C
76    SNESMonitorSolutionUpdate - Monitors progress of the SNES solvers by calling
77    VecView() for the UPDATE to the solution at each iteration.
78 
79    Collective on SNES
80 
81    Input Parameters:
82 +  snes - the SNES context
83 .  its - iteration number
84 .  fgnorm - 2-norm of residual
85 -  dummy - a viewer
86 
87    Options Database Keys:
88 .  -snes_monitor_solution_update [ascii binary draw][:filename][:viewer format] - plots update to solution at each iteration
89 
90    Level: intermediate
91 
92 .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
93 @*/
SNESMonitorSolutionUpdate(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)94 PetscErrorCode  SNESMonitorSolutionUpdate(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
95 {
96   PetscErrorCode ierr;
97   Vec            x;
98   PetscViewer    viewer = vf->viewer;
99 
100   PetscFunctionBegin;
101   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
102   ierr = SNESGetSolutionUpdate(snes,&x);CHKERRQ(ierr);
103   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
104   ierr = VecView(x,viewer);CHKERRQ(ierr);
105   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
106   PetscFunctionReturn(0);
107 }
108 
109 /*@C
110    KSPMonitorSNES - Print the residual norm of the nonlinear function at each iteration of the linear iterative solver.
111 
112    Collective on KSP
113 
114    Input Parameters:
115 +  ksp   - iterative context
116 .  n     - iteration number
117 .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
118 -  dummy - unused monitor context
119 
120    Level: intermediate
121 
122 .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
123 @*/
KSPMonitorSNES(KSP ksp,PetscInt n,PetscReal rnorm,void * dummy)124 PetscErrorCode  KSPMonitorSNES(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
125 {
126   PetscErrorCode ierr;
127   PetscViewer    viewer;
128   SNES           snes = (SNES) dummy;
129   Vec            snes_solution,work1,work2;
130   PetscReal      snorm;
131 
132   PetscFunctionBegin;
133   ierr = SNESGetSolution(snes,&snes_solution);CHKERRQ(ierr);
134   ierr = VecDuplicate(snes_solution,&work1);CHKERRQ(ierr);
135   ierr = VecDuplicate(snes_solution,&work2);CHKERRQ(ierr);
136   ierr = KSPBuildSolution(ksp,work1,NULL);CHKERRQ(ierr);
137   ierr = VecAYPX(work1,-1.0,snes_solution);CHKERRQ(ierr);
138   ierr = SNESComputeFunction(snes,work1,work2);CHKERRQ(ierr);
139   ierr = VecNorm(work2,NORM_2,&snorm);CHKERRQ(ierr);
140   ierr = VecDestroy(&work1);CHKERRQ(ierr);
141   ierr = VecDestroy(&work2);CHKERRQ(ierr);
142 
143   ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);CHKERRQ(ierr);
144   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
145   if (n == 0 && ((PetscObject)ksp)->prefix) {
146     ierr = PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);CHKERRQ(ierr);
147   }
148   ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Residual norm %5.3e KSP Residual norm %5.3e \n",n,(double)snorm,(double)rnorm);CHKERRQ(ierr);
149   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);CHKERRQ(ierr);
150   PetscFunctionReturn(0);
151 }
152 
153 #include <petscdraw.h>
154 
155 /*@C
156    KSPMonitorSNESLGResidualNormCreate - Creates a line graph context for use with
157    KSP to monitor convergence of preconditioned residual norms.
158 
159    Collective on KSP
160 
161    Input Parameters:
162 +  comm - communicator context
163 .  host - the X display to open, or null for the local machine
164 .  label - the title to put in the title bar
165 .  x, y - the screen coordinates of the upper left coordinate of
166           the window
167 -  m, n - the screen width and height in pixels
168 
169    Output Parameter:
170 .  draw - the drawing context
171 
172    Options Database Key:
173 .  -ksp_monitor_lg_residualnorm - Sets line graph monitor
174 
175    Notes:
176    Use KSPMonitorSNESLGResidualNormDestroy() to destroy this line graph; do not use PetscDrawLGDestroy().
177 
178    Level: intermediate
179 
180 .seealso: KSPMonitorSNESLGResidualNormDestroy(), KSPMonitorSet(), KSPMonitorSNESLGTrueResidualCreate()
181 @*/
KSPMonitorSNESLGResidualNormCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscObject ** objs)182 PetscErrorCode  KSPMonitorSNESLGResidualNormCreate(MPI_Comm comm,const char host[],const char label[],int x,int y,int m,int n,PetscObject **objs)
183 {
184   PetscDraw      draw;
185   PetscErrorCode ierr;
186   PetscDrawAxis  axis;
187   PetscDrawLG    lg;
188   const char     *names[] = {"Linear residual","Nonlinear residual"};
189 
190   PetscFunctionBegin;
191   ierr = PetscDrawCreate(comm,host,label,x,y,m,n,&draw);CHKERRQ(ierr);
192   ierr = PetscDrawSetFromOptions(draw);CHKERRQ(ierr);
193   ierr = PetscDrawLGCreate(draw,2,&lg);CHKERRQ(ierr);
194   ierr = PetscDrawLGSetLegend(lg,names);CHKERRQ(ierr);
195   ierr = PetscDrawLGSetFromOptions(lg);CHKERRQ(ierr);
196   ierr = PetscDrawLGGetAxis(lg,&axis);CHKERRQ(ierr);
197   ierr = PetscDrawAxisSetLabels(axis,"Convergence of Residual Norm","Iteration","Residual Norm");CHKERRQ(ierr);
198   ierr = PetscDrawDestroy(&draw);CHKERRQ(ierr);
199 
200   ierr = PetscMalloc1(2,objs);CHKERRQ(ierr);
201   (*objs)[1] = (PetscObject)lg;
202   PetscFunctionReturn(0);
203 }
204 
KSPMonitorSNESLGResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscObject * objs)205 PetscErrorCode  KSPMonitorSNESLGResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscObject *objs)
206 {
207   SNES           snes = (SNES) objs[0];
208   PetscDrawLG    lg   = (PetscDrawLG) objs[1];
209   PetscErrorCode ierr;
210   PetscReal      y[2];
211   Vec            snes_solution,work1,work2;
212 
213   PetscFunctionBegin;
214   if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
215   else y[0] = -15.0;
216 
217   ierr = SNESGetSolution(snes,&snes_solution);CHKERRQ(ierr);
218   ierr = VecDuplicate(snes_solution,&work1);CHKERRQ(ierr);
219   ierr = VecDuplicate(snes_solution,&work2);CHKERRQ(ierr);
220   ierr = KSPBuildSolution(ksp,work1,NULL);CHKERRQ(ierr);
221   ierr = VecAYPX(work1,-1.0,snes_solution);CHKERRQ(ierr);
222   ierr = SNESComputeFunction(snes,work1,work2);CHKERRQ(ierr);
223   ierr = VecNorm(work2,NORM_2,y+1);CHKERRQ(ierr);
224   if (y[1] > 0.0) y[1] = PetscLog10Real(y[1]);
225   else y[1] = -15.0;
226   ierr = VecDestroy(&work1);CHKERRQ(ierr);
227   ierr = VecDestroy(&work2);CHKERRQ(ierr);
228 
229   ierr = PetscDrawLGAddPoint(lg,NULL,y);CHKERRQ(ierr);
230   if (n < 20 || !(n % 5) || snes->reason) {
231     ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr);
232     ierr = PetscDrawLGSave(lg);CHKERRQ(ierr);
233   }
234   PetscFunctionReturn(0);
235 }
236 
237 /*@
238    KSPMonitorSNESLGResidualNormDestroy - Destroys a line graph context that was created
239    with KSPMonitorSNESLGResidualNormCreate().
240 
241    Collective on KSP
242 
243    Input Parameter:
244 .  draw - the drawing context
245 
246    Level: intermediate
247 
248 .seealso: KSPMonitorSNESLGResidualNormCreate(), KSPMonitorSNESLGTrueResidualDestroy(), KSPMonitorSet()
249 @*/
KSPMonitorSNESLGResidualNormDestroy(PetscObject ** objs)250 PetscErrorCode  KSPMonitorSNESLGResidualNormDestroy(PetscObject **objs)
251 {
252   PetscErrorCode ierr;
253   PetscDrawLG    lg = (PetscDrawLG) (*objs)[1];
254 
255   PetscFunctionBegin;
256   ierr = PetscDrawLGDestroy(&lg);CHKERRQ(ierr);
257   ierr = PetscFree(*objs);CHKERRQ(ierr);
258   PetscFunctionReturn(0);
259 }
260 
261 /*@C
262    SNESMonitorDefault - Monitors progress of the SNES solvers (default).
263 
264    Collective on SNES
265 
266    Input Parameters:
267 +  snes - the SNES context
268 .  its - iteration number
269 .  fgnorm - 2-norm of residual
270 -  vf - viewer and format structure
271 
272    Notes:
273    This routine prints the residual norm at each iteration.
274 
275    Level: intermediate
276 
277 .seealso: SNESMonitorSet(), SNESMonitorSolution()
278 @*/
SNESMonitorDefault(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)279 PetscErrorCode  SNESMonitorDefault(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
280 {
281   PetscErrorCode ierr;
282   PetscViewer    viewer = vf->viewer;
283 
284   PetscFunctionBegin;
285   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
286   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
287   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
288   ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %14.12e \n",its,(double)fgnorm);CHKERRQ(ierr);
289   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
290   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
291   PetscFunctionReturn(0);
292 }
293 
294 /*@C
295    SNESMonitorScaling - Monitors the largest value in each row of the Jacobian.
296 
297    Collective on SNES
298 
299    Input Parameters:
300 +  snes - the SNES context
301 .  its - iteration number
302 .  fgnorm - 2-norm of residual
303 -  vf - viewer and format structure
304 
305    Notes:
306    This routine prints the largest value in each row of the Jacobian
307 
308    Level: intermediate
309 
310 .seealso: SNESMonitorSet(), SNESMonitorSolution()
311 @*/
SNESMonitorScaling(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)312 PetscErrorCode  SNESMonitorScaling(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
313 {
314   PetscErrorCode ierr;
315   PetscViewer    viewer = vf->viewer;
316   KSP            ksp;
317   Mat            J;
318   Vec            v;
319 
320   PetscFunctionBegin;
321   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
322   ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
323   ierr = KSPGetOperators(ksp,&J,NULL);CHKERRQ(ierr);
324   ierr = MatCreateVecs(J,&v,NULL);CHKERRQ(ierr);
325   ierr = MatGetRowMaxAbs(J,v,NULL);CHKERRQ(ierr);
326   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
327   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
328   ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Jacobian maximum row entries \n");CHKERRQ(ierr);
329   ierr = VecView(v,viewer);CHKERRQ(ierr);
330   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
331   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
332   ierr = VecDestroy(&v);CHKERRQ(ierr);
333   PetscFunctionReturn(0);
334 }
335 
SNESMonitorJacUpdateSpectrum(SNES snes,PetscInt it,PetscReal fnorm,PetscViewerAndFormat * vf)336 PetscErrorCode SNESMonitorJacUpdateSpectrum(SNES snes,PetscInt it,PetscReal fnorm,PetscViewerAndFormat *vf)
337 {
338 #if defined(PETSC_HAVE_ESSL)
339   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - No support for ESSL Lapack Routines");
340 #else
341   Vec            X;
342   Mat            J,dJ,dJdense;
343   PetscErrorCode ierr;
344   PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);
345   PetscInt       n,i;
346   PetscBLASInt   nb = 0,lwork;
347   PetscReal      *eigr,*eigi;
348   PetscScalar    *work;
349   PetscScalar    *a;
350 
351   PetscFunctionBegin;
352   if (it == 0) PetscFunctionReturn(0);
353   /* create the difference between the current update and the current jacobian */
354   ierr = SNESGetSolution(snes,&X);CHKERRQ(ierr);
355   ierr = SNESGetJacobian(snes,NULL,&J,&func,NULL);CHKERRQ(ierr);
356   ierr = MatDuplicate(J,MAT_COPY_VALUES,&dJ);CHKERRQ(ierr);
357   ierr = SNESComputeJacobian(snes,X,dJ,dJ);CHKERRQ(ierr);
358   ierr = MatAXPY(dJ,-1.0,J,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
359 
360   /* compute the spectrum directly */
361   ierr  = MatConvert(dJ,MATSEQDENSE,MAT_INITIAL_MATRIX,&dJdense);CHKERRQ(ierr);
362   ierr  = MatGetSize(dJ,&n,NULL);CHKERRQ(ierr);
363   ierr  = PetscBLASIntCast(n,&nb);CHKERRQ(ierr);
364   lwork = 3*nb;
365   ierr  = PetscMalloc1(n,&eigr);CHKERRQ(ierr);
366   ierr  = PetscMalloc1(n,&eigi);CHKERRQ(ierr);
367   ierr  = PetscMalloc1(lwork,&work);CHKERRQ(ierr);
368   ierr  = MatDenseGetArray(dJdense,&a);CHKERRQ(ierr);
369 #if !defined(PETSC_USE_COMPLEX)
370   {
371     PetscBLASInt lierr;
372     ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
373     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,a,&nb,eigr,eigi,NULL,&nb,NULL,&nb,work,&lwork,&lierr));
374     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"geev() error %d",lierr);
375     ierr = PetscFPTrapPop();CHKERRQ(ierr);
376   }
377 #else
378   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
379 #endif
380   ierr = PetscPrintf(PetscObjectComm((PetscObject)snes),"Eigenvalues of J_%d - J_%d:\n",it,it-1);CHKERRQ(ierr);
381   for (i=0;i<n;i++) {
382     ierr = PetscPrintf(PetscObjectComm((PetscObject)snes),"%5d: %20.5g + %20.5gi\n",i,(double)eigr[i],(double)eigi[i]);CHKERRQ(ierr);
383   }
384   ierr = MatDenseRestoreArray(dJdense,&a);CHKERRQ(ierr);
385   ierr = MatDestroy(&dJ);CHKERRQ(ierr);
386   ierr = MatDestroy(&dJdense);CHKERRQ(ierr);
387   ierr = PetscFree(eigr);CHKERRQ(ierr);
388   ierr = PetscFree(eigi);CHKERRQ(ierr);
389   ierr = PetscFree(work);CHKERRQ(ierr);
390   PetscFunctionReturn(0);
391 #endif
392 }
393 
394 PETSC_INTERN PetscErrorCode  SNESMonitorRange_Private(SNES,PetscInt,PetscReal*);
395 
SNESMonitorRange_Private(SNES snes,PetscInt it,PetscReal * per)396 PetscErrorCode  SNESMonitorRange_Private(SNES snes,PetscInt it,PetscReal *per)
397 {
398   PetscErrorCode ierr;
399   Vec            resid;
400   PetscReal      rmax,pwork;
401   PetscInt       i,n,N;
402   PetscScalar    *r;
403 
404   PetscFunctionBegin;
405   ierr  = SNESGetFunction(snes,&resid,NULL,NULL);CHKERRQ(ierr);
406   ierr  = VecNorm(resid,NORM_INFINITY,&rmax);CHKERRQ(ierr);
407   ierr  = VecGetLocalSize(resid,&n);CHKERRQ(ierr);
408   ierr  = VecGetSize(resid,&N);CHKERRQ(ierr);
409   ierr  = VecGetArray(resid,&r);CHKERRQ(ierr);
410   pwork = 0.0;
411   for (i=0; i<n; i++) {
412     pwork += (PetscAbsScalar(r[i]) > .20*rmax);
413   }
414   ierr = MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
415   ierr = VecRestoreArray(resid,&r);CHKERRQ(ierr);
416   *per = *per/N;
417   PetscFunctionReturn(0);
418 }
419 
420 /*@C
421    SNESMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
422 
423    Collective on SNES
424 
425    Input Parameters:
426 +  snes   - iterative context
427 .  it    - iteration number
428 .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
429 -  dummy - unused monitor context
430 
431    Options Database Key:
432 .  -snes_monitor_range - Activates SNESMonitorRange()
433 
434    Level: intermediate
435 
436 .seealso: SNESMonitorSet(), SNESMonitorDefault(), SNESMonitorLGCreate()
437 @*/
SNESMonitorRange(SNES snes,PetscInt it,PetscReal rnorm,PetscViewerAndFormat * vf)438 PetscErrorCode  SNESMonitorRange(SNES snes,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *vf)
439 {
440   PetscErrorCode ierr;
441   PetscReal      perc,rel;
442   PetscViewer    viewer = vf->viewer;
443   /* should be in a MonitorRangeContext */
444   static PetscReal prev;
445 
446   PetscFunctionBegin;
447   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
448   if (!it) prev = rnorm;
449   ierr = SNESMonitorRange_Private(snes,it,&perc);CHKERRQ(ierr);
450 
451   rel  = (prev - rnorm)/prev;
452   prev = rnorm;
453   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
454   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
455   ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES preconditioned resid norm %14.12e Percent values above 20 percent of maximum %5.2f relative decrease %5.2e ratio %5.2e \n",it,(double)rnorm,(double)(100.0*perc),(double)rel,(double)(rel/perc));CHKERRQ(ierr);
456   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
457   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
458   PetscFunctionReturn(0);
459 }
460 
461 /*@C
462    SNESMonitorRatio - Monitors progress of the SNES solvers by printing the ratio
463    of residual norm at each iteration to the previous.
464 
465    Collective on SNES
466 
467    Input Parameters:
468 +  snes - the SNES context
469 .  its - iteration number
470 .  fgnorm - 2-norm of residual (or gradient)
471 -  dummy -  context of monitor
472 
473    Level: intermediate
474 
475    Notes:
476     Insure that SNESMonitorRatio() is called when you set this monitor
477 .seealso: SNESMonitorSet(), SNESMonitorSolution(), SNESMonitorRatio()
478 @*/
SNESMonitorRatio(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)479 PetscErrorCode  SNESMonitorRatio(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
480 {
481   PetscErrorCode          ierr;
482   PetscInt                len;
483   PetscReal               *history;
484   PetscViewer             viewer = vf->viewer;
485 
486   PetscFunctionBegin;
487   ierr = SNESGetConvergenceHistory(snes,&history,NULL,&len);CHKERRQ(ierr);
488   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
489   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
490   if (!its || !history || its > len) {
491     ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %14.12e \n",its,(double)fgnorm);CHKERRQ(ierr);
492   } else {
493     PetscReal ratio = fgnorm/history[its-1];
494     ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %14.12e %14.12e \n",its,(double)fgnorm,(double)ratio);CHKERRQ(ierr);
495   }
496   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
497   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
498   PetscFunctionReturn(0);
499 }
500 
501 /*@C
502    SNESMonitorRatioSetUp - Insures the SNES object is saving its history since this monitor needs access to it
503 
504    Collective on SNES
505 
506    Input Parameters:
507 +   snes - the SNES context
508 -   viewer - the PetscViewer object (ignored)
509 
510    Level: intermediate
511 
512 .seealso: SNESMonitorSet(), SNESMonitorSolution(), SNESMonitorDefault(), SNESMonitorRatio()
513 @*/
SNESMonitorRatioSetUp(SNES snes,PetscViewerAndFormat * vf)514 PetscErrorCode  SNESMonitorRatioSetUp(SNES snes,PetscViewerAndFormat *vf)
515 {
516   PetscErrorCode          ierr;
517   PetscReal               *history;
518 
519   PetscFunctionBegin;
520   ierr = SNESGetConvergenceHistory(snes,&history,NULL,NULL);CHKERRQ(ierr);
521   if (!history) {
522     ierr = SNESSetConvergenceHistory(snes,NULL,NULL,100,PETSC_TRUE);CHKERRQ(ierr);
523   }
524   PetscFunctionReturn(0);
525 }
526 
527 /* ---------------------------------------------------------------- */
528 /*
529      Default (short) SNES Monitor, same as SNESMonitorDefault() except
530   it prints fewer digits of the residual as the residual gets smaller.
531   This is because the later digits are meaningless and are often
532   different on different machines; by using this routine different
533   machines will usually generate the same output.
534 
535   Deprecated: Intentionally has no manual page
536 */
SNESMonitorDefaultShort(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)537 PetscErrorCode  SNESMonitorDefaultShort(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *vf)
538 {
539   PetscErrorCode ierr;
540   PetscViewer    viewer = vf->viewer;
541 
542   PetscFunctionBegin;
543   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
544   ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
545   ierr = PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
546   if (fgnorm > 1.e-9) {
547     ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %g \n",its,(double)fgnorm);CHKERRQ(ierr);
548   } else if (fgnorm > 1.e-11) {
549     ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %5.3e \n",its,(double)fgnorm);CHKERRQ(ierr);
550   } else {
551     ierr = PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm < 1.e-11\n",its);CHKERRQ(ierr);
552   }
553   ierr = PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
554   ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
555   PetscFunctionReturn(0);
556 }
557 
558 /*@C
559   SNESMonitorDefaultField - Monitors progress of the SNES solvers, separated into fields.
560 
561   Collective on SNES
562 
563   Input Parameters:
564 + snes   - the SNES context
565 . its    - iteration number
566 . fgnorm - 2-norm of residual
567 - ctx    - the PetscViewer
568 
569   Notes:
570   This routine uses the DM attached to the residual vector
571 
572   Level: intermediate
573 
574 .seealso: SNESMonitorSet(), SNESMonitorSolution(), SNESMonitorDefault()
575 @*/
SNESMonitorDefaultField(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)576 PetscErrorCode SNESMonitorDefaultField(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
577 {
578   PetscViewer    viewer = vf->viewer;
579   Vec            r;
580   DM             dm;
581   PetscReal      res[256];
582   PetscInt       tablevel;
583   PetscErrorCode ierr;
584 
585   PetscFunctionBegin;
586   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,4);
587   ierr = SNESGetFunction(snes, &r, NULL, NULL);CHKERRQ(ierr);
588   ierr = VecGetDM(r, &dm);CHKERRQ(ierr);
589   if (!dm) {ierr = SNESMonitorDefault(snes, its, fgnorm, vf);CHKERRQ(ierr);}
590   else {
591     PetscSection s, gs;
592     PetscInt     Nf, f;
593 
594     ierr = DMGetLocalSection(dm, &s);CHKERRQ(ierr);
595     ierr = DMGetGlobalSection(dm, &gs);CHKERRQ(ierr);
596     if (!s || !gs) {ierr = SNESMonitorDefault(snes, its, fgnorm, vf);CHKERRQ(ierr);}
597     ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr);
598     if (Nf > 256) SETERRQ1(PetscObjectComm((PetscObject) snes), PETSC_ERR_SUP, "Do not support %d fields > 256", Nf);
599     ierr = PetscSectionVecNorm(s, gs, r, NORM_2, res);CHKERRQ(ierr);
600     ierr = PetscObjectGetTabLevel((PetscObject) snes, &tablevel);CHKERRQ(ierr);
601     ierr = PetscViewerPushFormat(viewer,vf->format);CHKERRQ(ierr);
602     ierr = PetscViewerASCIIAddTab(viewer, tablevel);CHKERRQ(ierr);
603     ierr = PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);CHKERRQ(ierr);
604     for (f = 0; f < Nf; ++f) {
605       if (f) {ierr = PetscViewerASCIIPrintf(viewer, ", ");CHKERRQ(ierr);}
606       ierr = PetscViewerASCIIPrintf(viewer, "%14.12e", res[f]);CHKERRQ(ierr);
607     }
608     ierr = PetscViewerASCIIPrintf(viewer, "] \n");CHKERRQ(ierr);
609     ierr = PetscViewerASCIISubtractTab(viewer, tablevel);CHKERRQ(ierr);
610     ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
611   }
612   PetscFunctionReturn(0);
613 }
614 /* ---------------------------------------------------------------- */
615 /*@C
616    SNESConvergedDefault - Default onvergence test of the solvers for
617    systems of nonlinear equations.
618 
619    Collective on SNES
620 
621    Input Parameters:
622 +  snes - the SNES context
623 .  it - the iteration (0 indicates before any Newton steps)
624 .  xnorm - 2-norm of current iterate
625 .  snorm - 2-norm of current step
626 .  fnorm - 2-norm of function at current iterate
627 -  dummy - unused context
628 
629    Output Parameter:
630 .   reason  - one of
631 $  SNES_CONVERGED_FNORM_ABS       - (fnorm < abstol),
632 $  SNES_CONVERGED_SNORM_RELATIVE  - (snorm < stol*xnorm),
633 $  SNES_CONVERGED_FNORM_RELATIVE  - (fnorm < rtol*fnorm0),
634 $  SNES_DIVERGED_FUNCTION_COUNT   - (nfct > maxf),
635 $  SNES_DIVERGED_FNORM_NAN        - (fnorm == NaN),
636 $  SNES_CONVERGED_ITERATING       - (otherwise),
637 $  SNES_DIVERGED_DTOL             - (fnorm > divtol*snes->fnorm0)
638 
639    where
640 +    maxf - maximum number of function evaluations,  set with SNESSetTolerances()
641 .    nfct - number of function evaluations,
642 .    abstol - absolute function norm tolerance, set with SNESSetTolerances()
643 .    rtol - relative function norm tolerance, set with SNESSetTolerances()
644 .    divtol - divergence tolerance, set with SNESSetDivergenceTolerance()
645 -    fnorm0 - 2-norm of the function at the initial solution (initial guess; zeroth iteration)
646 
647   Options Database Keys:
648 +  -snes_stol - convergence tolerance in terms of the norm  of the change in the solution between steps
649 .  -snes_atol <abstol> - absolute tolerance of residual norm
650 .  -snes_rtol <rtol> - relative decrease in tolerance norm from the initial 2-norm of the solution
651 .  -snes_divergence_tolerance <divtol> - if the residual goes above divtol*rnorm0, exit with divergence
652 .  -snes_max_funcs <max_funcs> - maximum number of function evaluations
653 .  -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none
654 -  -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops
655 
656    Level: intermediate
657 
658 .seealso: SNESSetConvergenceTest(), SNESSetTolerances(), SNESSetDivergenceTolerance()
659 @*/
SNESConvergedDefault(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason * reason,void * dummy)660 PetscErrorCode  SNESConvergedDefault(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
661 {
662   PetscErrorCode ierr;
663 
664   PetscFunctionBegin;
665   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
666   PetscValidPointer(reason,6);
667 
668   *reason = SNES_CONVERGED_ITERATING;
669 
670   if (!it) {
671     /* set parameter for default relative tolerance convergence test */
672     snes->ttol   = fnorm*snes->rtol;
673     snes->rnorm0 = fnorm;
674   }
675   if (PetscIsInfOrNanReal(fnorm)) {
676     ierr    = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr);
677     *reason = SNES_DIVERGED_FNORM_NAN;
678   } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
679     ierr    = PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e\n",(double)fnorm,(double)snes->abstol);CHKERRQ(ierr);
680     *reason = SNES_CONVERGED_FNORM_ABS;
681   } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
682     ierr    = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr);
683     *reason = SNES_DIVERGED_FUNCTION_COUNT;
684   }
685 
686   if (it && !*reason) {
687     if (fnorm <= snes->ttol) {
688       ierr    = PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e (relative tolerance)\n",(double)fnorm,(double)snes->ttol);CHKERRQ(ierr);
689       *reason = SNES_CONVERGED_FNORM_RELATIVE;
690     } else if (snorm < snes->stol*xnorm) {
691       ierr    = PetscInfo3(snes,"Converged due to small update length: %14.12e < %14.12e * %14.12e\n",(double)snorm,(double)snes->stol,(double)xnorm);CHKERRQ(ierr);
692       *reason = SNES_CONVERGED_SNORM_RELATIVE;
693     } else if (snes->divtol > 0 && (fnorm > snes->divtol*snes->rnorm0)) {
694       ierr    = PetscInfo3(snes,"Diverged due to increase in function norm: %14.12e > %14.12e * %14.12e\n",(double)fnorm,(double)snes->divtol,(double)snes->rnorm0);CHKERRQ(ierr);
695       *reason = SNES_DIVERGED_DTOL;
696     }
697 
698   }
699   PetscFunctionReturn(0);
700 }
701 
702 /*@C
703    SNESConvergedSkip - Convergence test for SNES that NEVER returns as
704    converged, UNLESS the maximum number of iteration have been reached.
705 
706    Logically Collective on SNES
707 
708    Input Parameters:
709 +  snes - the SNES context
710 .  it - the iteration (0 indicates before any Newton steps)
711 .  xnorm - 2-norm of current iterate
712 .  snorm - 2-norm of current step
713 .  fnorm - 2-norm of function at current iterate
714 -  dummy - unused context
715 
716    Output Parameter:
717 .   reason  - SNES_CONVERGED_ITERATING, SNES_CONVERGED_ITS, or SNES_DIVERGED_FNORM_NAN
718 
719    Notes:
720    Convergence is then declared after a fixed number of iterations have been used.
721 
722    Level: advanced
723 
724 .seealso: SNESSetConvergenceTest()
725 @*/
SNESConvergedSkip(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason * reason,void * dummy)726 PetscErrorCode  SNESConvergedSkip(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
727 {
728   PetscErrorCode ierr;
729 
730   PetscFunctionBegin;
731   PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
732   PetscValidPointer(reason,6);
733 
734   *reason = SNES_CONVERGED_ITERATING;
735 
736   if (fnorm != fnorm) {
737     ierr    = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr);
738     *reason = SNES_DIVERGED_FNORM_NAN;
739   } else if (it == snes->max_its) {
740     *reason = SNES_CONVERGED_ITS;
741   }
742   PetscFunctionReturn(0);
743 }
744 
745 /*@C
746   SNESSetWorkVecs - Gets a number of work vectors.
747 
748   Input Parameters:
749 + snes  - the SNES context
750 - nw - number of work vectors to allocate
751 
752   Level: developer
753 
754 @*/
SNESSetWorkVecs(SNES snes,PetscInt nw)755 PetscErrorCode SNESSetWorkVecs(SNES snes,PetscInt nw)
756 {
757   DM             dm;
758   Vec            v;
759   PetscErrorCode ierr;
760 
761   PetscFunctionBegin;
762   if (snes->work) {ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr);}
763   snes->nwork = nw;
764 
765   ierr = SNESGetDM(snes, &dm);CHKERRQ(ierr);
766   ierr = DMGetGlobalVector(dm, &v);CHKERRQ(ierr);
767   ierr = VecDuplicateVecs(v,snes->nwork,&snes->work);CHKERRQ(ierr);
768   ierr = DMRestoreGlobalVector(dm, &v);CHKERRQ(ierr);
769   ierr = PetscLogObjectParents(snes,nw,snes->work);CHKERRQ(ierr);
770   PetscFunctionReturn(0);
771 }
772