Actual source code: ex2.c
petsc-3.9.4 2018-09-11
2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
3: This example employs a user-defined monitoring routine.\n\n";
5: /*T
6: Concepts: SNES^basic uniprocessor example
7: Concepts: SNES^setting a user-defined monitoring routine
8: Processors: 1
9: T*/
13: /*
14: Include "petscdraw.h" so that we can use PETSc drawing routines.
15: Include "petscsnes.h" so that we can use SNES solvers. Note that this
16: file automatically includes:
17: petscsys.h - base PETSc routines petscvec.h - vectors
18: petscmat.h - matrices
19: petscis.h - index sets petscksp.h - Krylov subspace methods
20: petscviewer.h - viewers petscpc.h - preconditioners
21: petscksp.h - linear solvers
22: */
24: #include <petscsnes.h>
26: /*
27: User-defined routines
28: */
29: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
30: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
31: extern PetscErrorCode FormInitialGuess(Vec);
32: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
34: /*
35: User-defined context for monitoring
36: */
37: typedef struct {
38: PetscViewer viewer;
39: } MonitorCtx;
41: int main(int argc,char **argv)
42: {
43: SNES snes; /* SNES context */
44: Vec x,r,F,U; /* vectors */
45: Mat J; /* Jacobian matrix */
46: MonitorCtx monP; /* monitoring context */
48: PetscInt its,n = 5,i,maxit,maxf;
49: PetscMPIInt size;
50: PetscScalar h,xp,v,none = -1.0;
51: PetscReal abstol,rtol,stol,norm;
53: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
54: MPI_Comm_size(PETSC_COMM_WORLD,&size);
55: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
56: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
57: h = 1.0/(n-1);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Create nonlinear solver context
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: SNESCreate(PETSC_COMM_WORLD,&snes);
65: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66: Create vector data structures; set function evaluation routine
67: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68: /*
69: Note that we form 1 vector from scratch and then duplicate as needed.
70: */
71: VecCreate(PETSC_COMM_WORLD,&x);
72: VecSetSizes(x,PETSC_DECIDE,n);
73: VecSetFromOptions(x);
74: VecDuplicate(x,&r);
75: VecDuplicate(x,&F);
76: VecDuplicate(x,&U);
78: /*
79: Set function evaluation routine and vector
80: */
81: SNESSetFunction(snes,r,FormFunction,(void*)F);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create matrix data structure; set Jacobian evaluation routine
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: MatCreate(PETSC_COMM_WORLD,&J);
89: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
90: MatSetFromOptions(J);
91: MatSeqAIJSetPreallocation(J,3,NULL);
93: /*
94: Set Jacobian matrix data structure and default Jacobian evaluation
95: routine. User can override with:
96: -snes_fd : default finite differencing approximation of Jacobian
97: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
98: (unless user explicitly sets preconditioner)
99: -snes_mf_operator : form preconditioning matrix as set by the user,
100: but use matrix-free approx for Jacobian-vector
101: products within Newton-Krylov method
102: */
104: SNESSetJacobian(snes,J,J,FormJacobian,NULL);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Customize nonlinear solver; set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Set an optional user-defined monitoring routine
112: */
113: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
114: SNESMonitorSet(snes,Monitor,&monP,0);
116: /*
117: Set names for some vectors to facilitate monitoring (optional)
118: */
119: PetscObjectSetName((PetscObject)x,"Approximate Solution");
120: PetscObjectSetName((PetscObject)U,"Exact Solution");
122: /*
123: Set SNES/KSP/KSP/PC runtime options, e.g.,
124: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
125: */
126: SNESSetFromOptions(snes);
128: /*
129: Print parameters used for convergence testing (optional) ... just
130: to demonstrate this routine; this information is also printed with
131: the option -snes_view
132: */
133: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
134: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Initialize application:
138: Store right-hand-side of PDE and exact solution
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: xp = 0.0;
142: for (i=0; i<n; i++) {
143: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
144: VecSetValues(F,1,&i,&v,INSERT_VALUES);
145: v = xp*xp*xp;
146: VecSetValues(U,1,&i,&v,INSERT_VALUES);
147: xp += h;
148: }
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Evaluate initial guess; then solve nonlinear system
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: /*
154: Note: The user should initialize the vector, x, with the initial guess
155: for the nonlinear solver prior to calling SNESSolve(). In particular,
156: to employ an initial guess of zero, the user should explicitly set
157: this vector to zero by calling VecSet().
158: */
159: FormInitialGuess(x);
160: SNESSolve(snes,NULL,x);
161: SNESGetIterationNumber(snes,&its);
162: PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Check solution and clean up
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: /*
169: Check the error
170: */
171: VecAXPY(x,none,U);
172: VecNorm(x,NORM_2,&norm);
173: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
176: /*
177: Free work space. All PETSc objects should be destroyed when they
178: are no longer needed.
179: */
180: VecDestroy(&x); VecDestroy(&r);
181: VecDestroy(&U); VecDestroy(&F);
182: MatDestroy(&J); SNESDestroy(&snes);
183: PetscViewerDestroy(&monP.viewer);
184: PetscFinalize();
185: return ierr;
186: }
187: /* ------------------------------------------------------------------- */
188: /*
189: FormInitialGuess - Computes initial guess.
191: Input/Output Parameter:
192: . x - the solution vector
193: */
194: PetscErrorCode FormInitialGuess(Vec x)
195: {
197: PetscScalar pfive = .50;
198: VecSet(x,pfive);
199: return 0;
200: }
201: /* ------------------------------------------------------------------- */
202: /*
203: FormFunction - Evaluates nonlinear function, F(x).
205: Input Parameters:
206: . snes - the SNES context
207: . x - input vector
208: . ctx - optional user-defined context, as set by SNESSetFunction()
210: Output Parameter:
211: . f - function vector
213: Note:
214: The user-defined context can contain any application-specific data
215: needed for the function evaluation (such as various parameters, work
216: vectors, and grid information). In this program the context is just
217: a vector containing the right-hand-side of the discretized PDE.
218: */
220: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
221: {
222: Vec g = (Vec)ctx;
223: const PetscScalar *xx,*gg;
224: PetscScalar *ff,d;
225: PetscErrorCode ierr;
226: PetscInt i,n;
228: /*
229: Get pointers to vector data.
230: - For default PETSc vectors, VecGetArray() returns a pointer to
231: the data array. Otherwise, the routine is implementation dependent.
232: - You MUST call VecRestoreArray() when you no longer need access to
233: the array.
234: */
235: VecGetArrayRead(x,&xx);
236: VecGetArray(f,&ff);
237: VecGetArrayRead(g,&gg);
239: /*
240: Compute function
241: */
242: VecGetSize(x,&n);
243: d = (PetscReal)(n - 1); d = d*d;
244: ff[0] = xx[0];
245: 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];
246: ff[n-1] = xx[n-1] - 1.0;
248: /*
249: Restore vectors
250: */
251: VecRestoreArrayRead(x,&xx);
252: VecRestoreArray(f,&ff);
253: VecRestoreArrayRead(g,&gg);
254: return 0;
255: }
256: /* ------------------------------------------------------------------- */
257: /*
258: FormJacobian - Evaluates Jacobian matrix.
260: Input Parameters:
261: . snes - the SNES context
262: . x - input vector
263: . dummy - optional user-defined context (not used here)
265: Output Parameters:
266: . jac - Jacobian matrix
267: . B - optionally different preconditioning matrix
269: */
271: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
272: {
273: const PetscScalar *xx;
274: PetscScalar A[3],d;
275: PetscErrorCode ierr;
276: PetscInt i,n,j[3];
278: /*
279: Get pointer to vector data
280: */
281: VecGetArrayRead(x,&xx);
283: /*
284: Compute Jacobian entries and insert into matrix.
285: - Note that in this case we set all elements for a particular
286: row at once.
287: */
288: VecGetSize(x,&n);
289: d = (PetscReal)(n - 1); d = d*d;
291: /*
292: Interior grid points
293: */
294: for (i=1; i<n-1; i++) {
295: j[0] = i - 1; j[1] = i; j[2] = i + 1;
296: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
297: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
298: }
300: /*
301: Boundary points
302: */
303: i = 0; A[0] = 1.0;
305: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
307: i = n-1; A[0] = 1.0;
309: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
311: /*
312: Restore vector
313: */
314: VecRestoreArrayRead(x,&xx);
316: /*
317: Assemble matrix
318: */
319: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
320: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
321: if (jac != B) {
322: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
323: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
324: }
325: return 0;
326: }
327: /* ------------------------------------------------------------------- */
328: /*
329: Monitor - User-defined monitoring routine that views the
330: current iterate with an x-window plot.
332: Input Parameters:
333: snes - the SNES context
334: its - iteration number
335: norm - 2-norm function value (may be estimated)
336: ctx - optional user-defined context for private data for the
337: monitor routine, as set by SNESMonitorSet()
339: Note:
340: See the manpage for PetscViewerDrawOpen() for useful runtime options,
341: such as -nox to deactivate all x-window output.
342: */
343: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
344: {
346: MonitorCtx *monP = (MonitorCtx*) ctx;
347: Vec x;
349: PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);
350: SNESGetSolution(snes,&x);
351: VecView(x,monP->viewer);
352: return 0;
353: }
356: /*TEST
358: test:
359: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
361: test:
362: suffix: 2
363: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
364: requires: !single
366: test:
367: suffix: 3
368: args: -nox -snes_monitor_cancel -snes_monitor_short -malloc no -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
370: TEST*/