Actual source code: ex1.c
petsc-3.3-p7 2013-05-11
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: typedef struct {
20: Vec xloc,rloc; /* local solution, residual vectors */
21: VecScatter scatter;
22: } AppCtx;
24: /*
25: User-defined routines
26: */
27: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
28: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
29: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
30: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
34: int main(int argc,char **argv)
35: {
36: SNES snes; /* nonlinear solver context */
37: KSP ksp; /* linear solver context */
38: PC pc; /* preconditioner context */
39: Vec x,r; /* solution, residual vectors */
40: Mat J; /* Jacobian matrix */
42: PetscInt its;
43: PetscMPIInt size,rank;
44: PetscScalar pfive = .5,*xx;
45: PetscBool flg;
46: AppCtx user; /* user-defined work context */
47: IS isglobal,islocal;
49: PetscInitialize(&argc,&argv,(char *)0,help);
50: MPI_Comm_size(PETSC_COMM_WORLD,&size);
51: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create nonlinear solver context
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56: SNESCreate(PETSC_COMM_WORLD,&snes);
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Create matrix and vector data structures; set corresponding routines
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: /*
62: Create vectors for solution and nonlinear function
63: */
64: VecCreate(PETSC_COMM_WORLD,&x);
65: VecSetSizes(x,PETSC_DECIDE,2);
66: VecSetFromOptions(x);
67: VecDuplicate(x,&r);
69: if (size > 1){
70: VecCreateSeq(PETSC_COMM_SELF,2,&user.xloc);
71: VecDuplicate(user.xloc,&user.rloc);
73: /* Create the scatter between the global x and local xloc */
74: ISCreateStride(MPI_COMM_SELF,2,0,1,&islocal);
75: ISCreateStride(MPI_COMM_SELF,2,0,1,&isglobal);
76: VecScatterCreate(x,isglobal,user.xloc,islocal,&user.scatter);
77: ISDestroy(&isglobal);
78: ISDestroy(&islocal);
79: }
81: /*
82: Create Jacobian matrix data structure
83: */
84: MatCreate(PETSC_COMM_WORLD,&J);
85: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
86: MatSetFromOptions(J);
87: MatSetUp(J);
89: PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
90: if (!flg) {
91: /*
92: Set function evaluation routine and vector.
93: */
94: SNESSetFunction(snes,r,FormFunction1,&user);
96: /*
97: Set Jacobian matrix data structure and Jacobian evaluation routine
98: */
99: SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
100: } else {
101: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This case is a uniprocessor example only!");
102: SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
103: SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
104: }
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Customize nonlinear solver; set runtime options
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: /*
110: Set linear solver defaults for this problem. By extracting the
111: KSP, KSP, and PC contexts from the SNES context, we can then
112: directly call any KSP, KSP, and PC routines to set various options.
113: */
114: SNESGetKSP(snes,&ksp);
115: KSPGetPC(ksp,&pc);
116: PCSetType(pc,PCNONE);
117: KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);
119: /*
120: Set SNES/KSP/KSP/PC runtime options, e.g.,
121: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
122: These options will override those specified above as long as
123: SNESSetFromOptions() is called _after_ any other customization
124: routines.
125: */
126: SNESSetFromOptions(snes);
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Evaluate initial guess; then solve nonlinear system
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: if (!flg) {
132: VecSet(x,pfive);
133: } else {
134: VecGetArray(x,&xx);
135: xx[0] = 2.0; xx[1] = 3.0;
136: VecRestoreArray(x,&xx);
137: }
138: /*
139: Note: The user should initialize the vector, x, with the initial guess
140: for the nonlinear solver prior to calling SNESSolve(). In particular,
141: to employ an initial guess of zero, the user should explicitly set
142: this vector to zero by calling VecSet().
143: */
145: SNESSolve(snes,PETSC_NULL,x);
146: SNESGetIterationNumber(snes,&its);
147: if (flg) {
148: Vec f;
149: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
150: SNESGetFunction(snes,&f,0,0);
151: VecView(r,PETSC_VIEWER_STDOUT_WORLD);
152: }
154: PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n",its);
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Free work space. All PETSc objects should be destroyed when they
158: are no longer needed.
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: VecDestroy(&x); VecDestroy(&r);
162: MatDestroy(&J); SNESDestroy(&snes);
163: if (size > 1){
164: VecDestroy(&user.xloc);
165: VecDestroy(&user.rloc);
166: VecScatterDestroy(&user.scatter);
167: }
168: PetscFinalize();
169: return 0;
170: }
171: /* ------------------------------------------------------------------- */
174: /*
175: FormFunction1 - Evaluates nonlinear function, F(x).
177: Input Parameters:
178: . snes - the SNES context
179: . x - input vector
180: . ctx - optional user-defined context
182: Output Parameter:
183: . f - function vector
184: */
185: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
186: {
187: PetscErrorCode ierr;
188: const PetscScalar *xx;
189: PetscScalar *ff;
190: AppCtx *user = (AppCtx*)ctx;
191: Vec xloc=user->xloc,floc=user->rloc;
192: VecScatter scatter=user->scatter;
193: MPI_Comm comm;
194: PetscMPIInt size,rank;
195: PetscInt rstart,rend;
197: PetscObjectGetComm((PetscObject)snes,&comm);
198: MPI_Comm_size(comm,&size);
199: MPI_Comm_rank(comm,&rank);
200: if (size > 1){
201: /*
202: This is a ridiculous case for testing intermidiate steps from sequential
203: code development to parallel implementation.
204: (1) scatter x into a sequetial vector;
205: (2) each process evaluates all values of floc;
206: (3) scatter floc back to the parallel f.
207: */
208: VecScatterBegin(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
209: VecScatterEnd(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
211: VecGetOwnershipRange(f,&rstart,&rend);
212: VecGetArrayRead(xloc,&xx);
213: VecGetArray(floc,&ff);
214: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
215: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
216: VecRestoreArray(floc,&ff);
217: VecRestoreArrayRead(xloc,&xx);
219: VecScatterBegin(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
220: VecScatterEnd(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
221: } else {
222: /*
223: Get pointers to vector data.
224: - For default PETSc vectors, VecGetArray() returns a pointer to
225: the data array. Otherwise, the routine is implementation dependent.
226: - You MUST call VecRestoreArray() when you no longer need access to
227: the array.
228: */
229: VecGetArrayRead(x,&xx);
230: VecGetArray(f,&ff);
232: /* Compute function */
233: ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
234: ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
236: /* Restore vectors */
237: VecRestoreArrayRead(x,&xx);
238: VecRestoreArray(f,&ff);
239: }
240: return 0;
241: }
242: /* ------------------------------------------------------------------- */
245: /*
246: FormJacobian1 - Evaluates Jacobian matrix.
248: Input Parameters:
249: . snes - the SNES context
250: . x - input vector
251: . dummy - optional user-defined context (not used here)
253: Output Parameters:
254: . jac - Jacobian matrix
255: . B - optionally different preconditioning matrix
256: . flag - flag indicating matrix structure
257: */
258: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
259: {
260: PetscScalar *xx,A[4];
262: PetscInt idx[2] = {0,1};
264: /*
265: Get pointer to vector data
266: */
267: VecGetArray(x,&xx);
269: /*
270: Compute Jacobian entries and insert into matrix.
271: - Since this is such a small problem, we set all entries for
272: the matrix at once.
273: */
274: A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
275: A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
276: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
277: *flag = SAME_NONZERO_PATTERN;
279: /*
280: Restore vector
281: */
282: VecRestoreArray(x,&xx);
284: /*
285: Assemble matrix
286: */
287: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
288: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
289: if (*jac != *B){
290: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
291: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
292: }
293: return 0;
294: }
296: /* ------------------------------------------------------------------- */
299: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
300: {
302: PetscScalar *xx,*ff;
304: /*
305: Get pointers to vector data.
306: - For default PETSc vectors, VecGetArray() returns a pointer to
307: the data array. Otherwise, the routine is implementation dependent.
308: - You MUST call VecRestoreArray() when you no longer need access to
309: the array.
310: */
311: VecGetArray(x,&xx);
312: VecGetArray(f,&ff);
314: /*
315: Compute function
316: */
317: ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
318: ff[1] = xx[1];
320: /*
321: Restore vectors
322: */
323: VecRestoreArray(x,&xx);
324: VecRestoreArray(f,&ff);
325: return 0;
326: }
327: /* ------------------------------------------------------------------- */
330: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
331: {
332: PetscScalar *xx,A[4];
334: PetscInt idx[2] = {0,1};
336: /*
337: Get pointer to vector data
338: */
339: VecGetArray(x,&xx);
341: /*
342: Compute Jacobian entries and insert into matrix.
343: - Since this is such a small problem, we set all entries for
344: the matrix at once.
345: */
346: A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
347: A[2] = 0.0; A[3] = 1.0;
348: MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
349: *flag = SAME_NONZERO_PATTERN;
351: /*
352: Restore vector
353: */
354: VecRestoreArray(x,&xx);
356: /*
357: Assemble matrix
358: */
359: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
360: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
361: if (*jac != *B){
362: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
363: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
364: }
365: return 0;
366: }