Actual source code: ex1.c
2: static char help[] = "Newton's method for a two-variable system, sequential.\n\n";
4: /*T
5: Concepts: SNES^basic example
6: T*/
10: /*
11: Include "petscsnes.h" so that we can use SNES solvers. Note that this
12: file automatically includes:
13: petscsys.h - base PETSc routines petscvec.h - vectors
14: petscmat.h - matrices
15: petscis.h - index sets petscksp.h - Krylov subspace methods
16: petscviewer.h - viewers petscpc.h - preconditioners
17: petscksp.h - linear solvers
18: */
This examples solves either
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6}
\end{equation}
or if the {\tt -hard} options is given
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
\end{equation}
29: #include <petscsnes.h>
31: /*
32: User-defined routines
33: */
34: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
35: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
36: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
37: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
39: int main(int argc,char **argv)
40: {
41: SNES snes; /* nonlinear solver context */
42: KSP ksp; /* linear solver context */
43: PC pc; /* preconditioner context */
44: Vec x,r; /* solution, residual vectors */
45: Mat J; /* Jacobian matrix */
47: PetscMPIInt size;
48: PetscScalar pfive = .5,*xx;
49: PetscBool flg;
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_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create nonlinear solver context
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: SNESCreate(PETSC_COMM_WORLD,&snes);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create matrix and vector data structures; set corresponding routines
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: /*
64: Create vectors for solution and nonlinear function
65: */
66: VecCreate(PETSC_COMM_WORLD,&x);
67: VecSetSizes(x,PETSC_DECIDE,2);
68: VecSetFromOptions(x);
69: VecDuplicate(x,&r);
71: /*
72: Create Jacobian matrix data structure
73: */
74: MatCreate(PETSC_COMM_WORLD,&J);
75: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
76: MatSetFromOptions(J);
77: MatSetUp(J);
79: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
80: if (!flg) {
81: /*
82: Set function evaluation routine and vector.
83: */
84: SNESSetFunction(snes,r,FormFunction1,NULL);
86: /*
87: Set Jacobian matrix data structure and Jacobian evaluation routine
88: */
89: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
90: } else {
91: SNESSetFunction(snes,r,FormFunction2,NULL);
92: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
93: }
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Customize nonlinear solver; set runtime options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: /*
99: Set linear solver defaults for this problem. By extracting the
100: KSP and PC contexts from the SNES context, we can then
101: directly call any KSP and PC routines to set various options.
102: */
103: SNESGetKSP(snes,&ksp);
104: KSPGetPC(ksp,&pc);
105: PCSetType(pc,PCNONE);
106: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
108: /*
109: Set SNES/KSP/KSP/PC runtime options, e.g.,
110: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
111: These options will override those specified above as long as
112: SNESSetFromOptions() is called _after_ any other customization
113: routines.
114: */
115: SNESSetFromOptions(snes);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Evaluate initial guess; then solve nonlinear system
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: if (!flg) {
121: VecSet(x,pfive);
122: } else {
123: VecGetArray(x,&xx);
124: xx[0] = 2.0; xx[1] = 3.0;
125: VecRestoreArray(x,&xx);
126: }
127: /*
128: Note: The user should initialize the vector, x, with the initial guess
129: for the nonlinear solver prior to calling SNESSolve(). In particular,
130: to employ an initial guess of zero, the user should explicitly set
131: this vector to zero by calling VecSet().
132: */
134: SNESSolve(snes,NULL,x);
135: if (flg) {
136: Vec f;
137: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
138: SNESGetFunction(snes,&f,0,0);
139: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
140: }
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Free work space. All PETSc objects should be destroyed when they
145: are no longer needed.
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: VecDestroy(&x); VecDestroy(&r);
149: MatDestroy(&J); SNESDestroy(&snes);
150: PetscFinalize();
151: return ierr;
152: }
153: /* ------------------------------------------------------------------- */
154: /*
155: FormFunction1 - Evaluates nonlinear function, F(x).
157: Input Parameters:
158: . snes - the SNES context
159: . x - input vector
160: . ctx - optional user-defined context
162: Output Parameter:
163: . f - function vector
164: */
165: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
166: {
167: PetscErrorCode ierr;
168: const PetscScalar *xx;
169: PetscScalar *ff;
171: /*
172: Get pointers to vector data.
173: - For default PETSc vectors, VecGetArray() returns a pointer to
174: the data array. Otherwise, the routine is implementation dependent.
175: - You MUST call VecRestoreArray() when you no longer need access to
176: the array.
177: */
178: VecGetArrayRead(x,&xx);
179: VecGetArray(f,&ff);
181: /* Compute function */
182: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
183: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
185: /* Restore vectors */
186: VecRestoreArrayRead(x,&xx);
187: VecRestoreArray(f,&ff);
188: return 0;
189: }
190: /* ------------------------------------------------------------------- */
191: /*
192: FormJacobian1 - Evaluates Jacobian matrix.
194: Input Parameters:
195: . snes - the SNES context
196: . x - input vector
197: . dummy - optional user-defined context (not used here)
199: Output Parameters:
200: . jac - Jacobian matrix
201: . B - optionally different preconditioning matrix
202: . flag - flag indicating matrix structure
203: */
204: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
205: {
206: const PetscScalar *xx;
207: PetscScalar A[4];
208: PetscErrorCode ierr;
209: PetscInt idx[2] = {0,1};
211: /*
212: Get pointer to vector data
213: */
214: VecGetArrayRead(x,&xx);
216: /*
217: Compute Jacobian entries and insert into matrix.
218: - Since this is such a small problem, we set all entries for
219: the matrix at once.
220: */
221: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
222: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
223: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
225: /*
226: Restore vector
227: */
228: VecRestoreArrayRead(x,&xx);
230: /*
231: Assemble matrix
232: */
233: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
234: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
235: if (jac != B) {
236: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
237: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
238: }
239: return 0;
240: }
242: /* ------------------------------------------------------------------- */
243: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
244: {
245: PetscErrorCode ierr;
246: const PetscScalar *xx;
247: PetscScalar *ff;
249: /*
250: Get pointers to vector data.
251: - For default PETSc vectors, VecGetArray() returns a pointer to
252: the data array. Otherwise, the routine is implementation dependent.
253: - You MUST call VecRestoreArray() when you no longer need access to
254: the array.
255: */
256: VecGetArrayRead(x,&xx);
257: VecGetArray(f,&ff);
259: /*
260: Compute function
261: */
262: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
263: ff[1] = xx[1];
265: /*
266: Restore vectors
267: */
268: VecRestoreArrayRead(x,&xx);
269: VecRestoreArray(f,&ff);
270: return 0;
271: }
272: /* ------------------------------------------------------------------- */
273: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
274: {
275: const PetscScalar *xx;
276: PetscScalar A[4];
277: PetscErrorCode ierr;
278: PetscInt idx[2] = {0,1};
280: /*
281: Get pointer to vector data
282: */
283: VecGetArrayRead(x,&xx);
285: /*
286: Compute Jacobian entries and insert into matrix.
287: - Since this is such a small problem, we set all entries for
288: the matrix at once.
289: */
290: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
291: A[2] = 0.0; A[3] = 1.0;
292: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
294: /*
295: Restore vector
296: */
297: VecRestoreArrayRead(x,&xx);
299: /*
300: Assemble matrix
301: */
302: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
303: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
304: if (jac != B) {
305: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
306: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
307: }
308: return 0;
309: }
314: /*TEST
316: test:
317: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
318: requires: !single
320: test:
321: suffix: 2
322: requires: !single
323: args: -snes_monitor_short
324: output_file: output/ex1_1.out
326: test:
327: suffix: 3
328: args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short
329: requires: !single
330: output_file: output/ex1_1.out
332: test:
333: suffix: 4
334: args: -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short
335: requires: !single
336: output_file: output/ex1_1.out
338: test:
339: suffix: 5
340: args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short
341: requires: !single
342: output_file: output/ex1_1.out
344: test:
345: suffix: 6
346: args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short
347: requires: !single
348: output_file: output/ex1_1.out
350: test:
351: suffix: X
352: args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
353: requires: !single x
355: TEST*/