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*/
8: /*
9: Include "petscsnes.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petscsys.h - base PETSc routines petscvec.h - vectors
12: petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
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}
27: #include <petscsnes.h>
29: /*
30: User-defined routines
31: */
32: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
33: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
34: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
35: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
37: int main(int argc,char **argv)
38: {
39: SNES snes; /* nonlinear solver context */
40: KSP ksp; /* linear solver context */
41: PC pc; /* preconditioner context */
42: Vec x,r; /* solution, residual vectors */
43: Mat J; /* Jacobian matrix */
44: PetscMPIInt size;
45: PetscScalar pfive = .5,*xx;
46: PetscBool flg;
48: PetscInitialize(&argc,&argv,(char*)0,help);
49: MPI_Comm_size(PETSC_COMM_WORLD,&size);
52: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: Create nonlinear solver context
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: SNESCreate(PETSC_COMM_WORLD,&snes);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create matrix and vector data structures; set corresponding routines
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: /*
61: Create vectors for solution and nonlinear function
62: */
63: VecCreate(PETSC_COMM_WORLD,&x);
64: VecSetSizes(x,PETSC_DECIDE,2);
65: VecSetFromOptions(x);
66: VecDuplicate(x,&r);
68: /*
69: Create Jacobian matrix data structure
70: */
71: MatCreate(PETSC_COMM_WORLD,&J);
72: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
73: MatSetFromOptions(J);
74: MatSetUp(J);
76: PetscOptionsHasName(NULL,NULL,"-hard",&flg);
77: if (!flg) {
78: /*
79: Set function evaluation routine and vector.
80: */
81: SNESSetFunction(snes,r,FormFunction1,NULL);
83: /*
84: Set Jacobian matrix data structure and Jacobian evaluation routine
85: */
86: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
87: } else {
88: SNESSetFunction(snes,r,FormFunction2,NULL);
89: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
90: }
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Customize nonlinear solver; set runtime options
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: /*
96: Set linear solver defaults for this problem. By extracting the
97: KSP and PC contexts from the SNES context, we can then
98: directly call any KSP and PC routines to set various options.
99: */
100: SNESGetKSP(snes,&ksp);
101: KSPGetPC(ksp,&pc);
102: PCSetType(pc,PCNONE);
103: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
105: /*
106: Set SNES/KSP/KSP/PC runtime options, e.g.,
107: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
108: These options will override those specified above as long as
109: SNESSetFromOptions() is called _after_ any other customization
110: routines.
111: */
112: SNESSetFromOptions(snes);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Evaluate initial guess; then solve nonlinear system
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: if (!flg) {
118: VecSet(x,pfive);
119: } else {
120: VecGetArray(x,&xx);
121: xx[0] = 2.0; xx[1] = 3.0;
122: VecRestoreArray(x,&xx);
123: }
124: /*
125: Note: The user should initialize the vector, x, with the initial guess
126: for the nonlinear solver prior to calling SNESSolve(). In particular,
127: to employ an initial guess of zero, the user should explicitly set
128: this vector to zero by calling VecSet().
129: */
131: SNESSolve(snes,NULL,x);
132: if (flg) {
133: Vec f;
134: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
135: SNESGetFunction(snes,&f,0,0);
136: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
137: }
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Free work space. All PETSc objects should be destroyed when they
141: are no longer needed.
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: VecDestroy(&x)); PetscCall(VecDestroy(&r);
145: MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
146: PetscFinalize();
147: return 0;
148: }
149: /* ------------------------------------------------------------------- */
150: /*
151: FormFunction1 - Evaluates nonlinear function, F(x).
153: Input Parameters:
154: . snes - the SNES context
155: . x - input vector
156: . ctx - optional user-defined context
158: Output Parameter:
159: . f - function vector
160: */
161: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
162: {
163: const PetscScalar *xx;
164: PetscScalar *ff;
166: /*
167: Get pointers to vector data.
168: - For default PETSc vectors, VecGetArray() returns a pointer to
169: the data array. Otherwise, the routine is implementation dependent.
170: - You MUST call VecRestoreArray() when you no longer need access to
171: the array.
172: */
173: VecGetArrayRead(x,&xx);
174: VecGetArray(f,&ff);
176: /* Compute function */
177: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
178: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
180: /* Restore vectors */
181: VecRestoreArrayRead(x,&xx);
182: VecRestoreArray(f,&ff);
183: return 0;
184: }
185: /* ------------------------------------------------------------------- */
186: /*
187: FormJacobian1 - Evaluates Jacobian matrix.
189: Input Parameters:
190: . snes - the SNES context
191: . x - input vector
192: . dummy - optional user-defined context (not used here)
194: Output Parameters:
195: . jac - Jacobian matrix
196: . B - optionally different preconditioning matrix
197: . flag - flag indicating matrix structure
198: */
199: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
200: {
201: const PetscScalar *xx;
202: PetscScalar A[4];
203: PetscInt idx[2] = {0,1};
205: /*
206: Get pointer to vector data
207: */
208: VecGetArrayRead(x,&xx);
210: /*
211: Compute Jacobian entries and insert into matrix.
212: - Since this is such a small problem, we set all entries for
213: the matrix at once.
214: */
215: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
216: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
217: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
219: /*
220: Restore vector
221: */
222: VecRestoreArrayRead(x,&xx);
224: /*
225: Assemble matrix
226: */
227: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
229: if (jac != B) {
230: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
231: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
232: }
233: return 0;
234: }
236: /* ------------------------------------------------------------------- */
237: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
238: {
239: const PetscScalar *xx;
240: PetscScalar *ff;
242: /*
243: Get pointers to vector data.
244: - For default PETSc vectors, VecGetArray() returns a pointer to
245: the data array. Otherwise, the routine is implementation dependent.
246: - You MUST call VecRestoreArray() when you no longer need access to
247: the array.
248: */
249: VecGetArrayRead(x,&xx);
250: VecGetArray(f,&ff);
252: /*
253: Compute function
254: */
255: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
256: ff[1] = xx[1];
258: /*
259: Restore vectors
260: */
261: VecRestoreArrayRead(x,&xx);
262: VecRestoreArray(f,&ff);
263: return 0;
264: }
265: /* ------------------------------------------------------------------- */
266: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
267: {
268: const PetscScalar *xx;
269: PetscScalar A[4];
270: PetscInt idx[2] = {0,1};
272: /*
273: Get pointer to vector data
274: */
275: VecGetArrayRead(x,&xx);
277: /*
278: Compute Jacobian entries and insert into matrix.
279: - Since this is such a small problem, we set all entries for
280: the matrix at once.
281: */
282: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
283: A[2] = 0.0; A[3] = 1.0;
284: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
286: /*
287: Restore vector
288: */
289: VecRestoreArrayRead(x,&xx);
291: /*
292: Assemble matrix
293: */
294: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
296: if (jac != B) {
297: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
298: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
299: }
300: return 0;
301: }
303: /*TEST
305: test:
306: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
307: requires: !single
309: test:
310: suffix: 2
311: requires: !single
312: args: -snes_monitor_short
313: output_file: output/ex1_1.out
315: test:
316: suffix: 3
317: args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab -snes_monitor_short
318: requires: !single
319: output_file: output/ex1_1.out
321: test:
322: suffix: 4
323: args: -ksp_view_solution ascii:ex1_2_sol.tmp::append -snes_monitor_short
324: requires: !single
325: output_file: output/ex1_1.out
327: test:
328: suffix: 5
329: args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append -snes_monitor_short
330: requires: !single
331: output_file: output/ex1_1.out
333: test:
334: suffix: 6
335: args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append -snes_monitor_short
336: requires: !single
337: output_file: output/ex1_1.out
339: test:
340: suffix: X
341: args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
342: requires: !single x
344: TEST*/