Actual source code: ex1.c
petsc-3.4.5 2014-06-29
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: */
17: #include <petscsnes.h>
19: /*
20: User-defined routines
21: */
22: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
23: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
24: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
25: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
29: int main(int argc,char **argv)
30: {
31: SNES snes; /* nonlinear solver context */
32: KSP ksp; /* linear solver context */
33: PC pc; /* preconditioner context */
34: Vec x,r; /* solution, residual vectors */
35: Mat J; /* Jacobian matrix */
37: PetscInt its;
38: PetscMPIInt size;
39: PetscScalar pfive = .5,*xx;
40: PetscBool flg;
42: PetscInitialize(&argc,&argv,(char*)0,help);
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Create nonlinear solver context
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: SNESCreate(PETSC_COMM_WORLD,&snes);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create matrix and vector data structures; set corresponding routines
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: /*
55: Create vectors for solution and nonlinear function
56: */
57: VecCreate(PETSC_COMM_WORLD,&x);
58: VecSetSizes(x,PETSC_DECIDE,2);
59: VecSetFromOptions(x);
60: VecDuplicate(x,&r);
62: /*
63: Create Jacobian matrix data structure
64: */
65: MatCreate(PETSC_COMM_WORLD,&J);
66: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
67: MatSetFromOptions(J);
68: MatSetUp(J);
70: PetscOptionsHasName(NULL,"-hard",&flg);
71: if (!flg) {
72: /*
73: Set function evaluation routine and vector.
74: */
75: SNESSetFunction(snes,r,FormFunction1,NULL);
77: /*
78: Set Jacobian matrix data structure and Jacobian evaluation routine
79: */
80: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
81: } else {
82: SNESSetFunction(snes,r,FormFunction2,NULL);
83: SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
84: }
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Customize nonlinear solver; set runtime options
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: /*
90: Set linear solver defaults for this problem. By extracting the
91: KSP and PC contexts from the SNES context, we can then
92: directly call any KSP and PC routines to set various options.
93: */
94: SNESGetKSP(snes,&ksp);
95: KSPGetPC(ksp,&pc);
96: PCSetType(pc,PCNONE);
97: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
99: /*
100: Set SNES/KSP/KSP/PC runtime options, e.g.,
101: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
102: These options will override those specified above as long as
103: SNESSetFromOptions() is called _after_ any other customization
104: routines.
105: */
106: SNESSetFromOptions(snes);
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Evaluate initial guess; then solve nonlinear system
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: if (!flg) {
112: VecSet(x,pfive);
113: } else {
114: VecGetArray(x,&xx);
115: xx[0] = 2.0; xx[1] = 3.0;
116: VecRestoreArray(x,&xx);
117: }
118: /*
119: Note: The user should initialize the vector, x, with the initial guess
120: for the nonlinear solver prior to calling SNESSolve(). In particular,
121: to employ an initial guess of zero, the user should explicitly set
122: this vector to zero by calling VecSet().
123: */
125: SNESSolve(snes,NULL,x);
126: SNESGetIterationNumber(snes,&its);
127: if (flg) {
128: Vec f;
129: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
130: SNESGetFunction(snes,&f,0,0);
131: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
132: }
134: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Free work space. All PETSc objects should be destroyed when they
138: are no longer needed.
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: VecDestroy(&x); VecDestroy(&r);
142: MatDestroy(&J); SNESDestroy(&snes);
143: PetscFinalize();
144: return 0;
145: }
146: /* ------------------------------------------------------------------- */
149: /*
150: FormFunction1 - Evaluates nonlinear function, F(x).
152: Input Parameters:
153: . snes - the SNES context
154: . x - input vector
155: . ctx - optional user-defined context
157: Output Parameter:
158: . f - function vector
159: */
160: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
161: {
162: PetscErrorCode ierr;
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: /* ------------------------------------------------------------------- */
188: /*
189: FormJacobian1 - Evaluates Jacobian matrix.
191: Input Parameters:
192: . snes - the SNES context
193: . x - input vector
194: . dummy - optional user-defined context (not used here)
196: Output Parameters:
197: . jac - Jacobian matrix
198: . B - optionally different preconditioning matrix
199: . flag - flag indicating matrix structure
200: */
201: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
202: {
203: const PetscScalar *xx;
204: PetscScalar A[4];
205: PetscErrorCode ierr;
206: PetscInt idx[2] = {0,1};
208: /*
209: Get pointer to vector data
210: */
211: VecGetArrayRead(x,&xx);
213: /*
214: Compute Jacobian entries and insert into matrix.
215: - Since this is such a small problem, we set all entries for
216: the matrix at once.
217: */
218: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
219: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
220: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
221: *flag = SAME_NONZERO_PATTERN;
223: /*
224: Restore vector
225: */
226: VecRestoreArrayRead(x,&xx);
228: /*
229: Assemble matrix
230: */
231: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
232: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
233: if (*jac != *B) {
234: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
235: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
236: }
237: return 0;
238: }
240: /* ------------------------------------------------------------------- */
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: /* ------------------------------------------------------------------- */
275: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
276: {
277: const PetscScalar *xx;
278: PetscScalar A[4];
279: PetscErrorCode ierr;
280: PetscInt idx[2] = {0,1};
282: /*
283: Get pointer to vector data
284: */
285: VecGetArrayRead(x,&xx);
287: /*
288: Compute Jacobian entries and insert into matrix.
289: - Since this is such a small problem, we set all entries for
290: the matrix at once.
291: */
292: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
293: A[2] = 0.0; A[3] = 1.0;
294: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
295: *flag = SAME_NONZERO_PATTERN;
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: }