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