Actual source code: ex1.c
petsc-3.6.1 2015-08-06
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,void*);
23: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
24: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,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,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);
222: /*
223: Restore vector
224: */
225: VecRestoreArrayRead(x,&xx);
227: /*
228: Assemble matrix
229: */
230: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
231: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
232: if (jac != B) {
233: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
234: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
235: }
236: return 0;
237: }
239: /* ------------------------------------------------------------------- */
242: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
243: {
244: PetscErrorCode ierr;
245: const PetscScalar *xx;
246: PetscScalar *ff;
248: /*
249: Get pointers to vector data.
250: - For default PETSc vectors, VecGetArray() returns a pointer to
251: the data array. Otherwise, the routine is implementation dependent.
252: - You MUST call VecRestoreArray() when you no longer need access to
253: the array.
254: */
255: VecGetArrayRead(x,&xx);
256: VecGetArray(f,&ff);
258: /*
259: Compute function
260: */
261: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
262: ff[1] = xx[1];
264: /*
265: Restore vectors
266: */
267: VecRestoreArrayRead(x,&xx);
268: VecRestoreArray(f,&ff);
269: return 0;
270: }
271: /* ------------------------------------------------------------------- */
274: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
275: {
276: const PetscScalar *xx;
277: PetscScalar A[4];
278: PetscErrorCode ierr;
279: PetscInt idx[2] = {0,1};
281: /*
282: Get pointer to vector data
283: */
284: VecGetArrayRead(x,&xx);
286: /*
287: Compute Jacobian entries and insert into matrix.
288: - Since this is such a small problem, we set all entries for
289: the matrix at once.
290: */
291: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
292: A[2] = 0.0; A[3] = 1.0;
293: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
295: /*
296: Restore vector
297: */
298: VecRestoreArrayRead(x,&xx);
300: /*
301: Assemble matrix
302: */
303: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
304: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
305: if (jac != B) {
306: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
307: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
308: }
309: return 0;
310: }