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