Actual source code: ex1.c
petsc-3.11.4 2019-09-28
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*/