Actual source code: ex45.c
petsc-3.3-p7 2013-05-11
2: static char help[] = "u`` + u^{2} = f. \n\
3: Demonstrates the use of Newton-Krylov methods \n\
4: with different jacobian evaluation routines and matrix colorings. \n\
5: Modified from ex6.c \n\
6: Input arguments are:\n\
7: -snes_jacobian_default : Jacobian using finite differences. Slow and expensive, not take advantage of sparsity \n\
8: -fd_jacobian_coloring: Jacobian using finite differences with matrix coloring\n\
9: -my_jacobian_struct: use user-provided Jacobian data structure to create matcoloring context \n\n";
11: /*
12: Example:
13: ./ex45 -n 10 -snes_monitor -ksp_monitor
14: ./ex45 -n 10 -snes_monitor -ksp_monitor -snes_jacobian_default -pc_type jacobi
15: ./ex45 -n 10 -snes_monitor -ksp_monitor -snes_jacobian_default -pc_type ilu
16: ./ex45 -n 10 -snes_jacobian_default -log_summary |grep SNESFunctionEval
17: ./ex45 -n 10 -snes_jacobian_default -fd_jacobian_coloring -my_jacobian_struct -log_summary |grep SNESFunctionEval
18: */
20: #include <petscsnes.h>
22: /*
23: User-defined routines
24: */
25: PetscErrorCode MyJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
26: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
27: PetscErrorCode MyApproxJacobianStructure(Mat *,void *);
29: int main(int argc,char **argv)
30: {
31: SNES snes; /* SNES context */
32: Vec x,r,F; /* vectors */
33: Mat J,JPrec; /* Jacobian,preconditioner matrices */
35: PetscInt it,n = 5,i;
36: PetscMPIInt size;
37: PetscReal h,xp = 0.0;
38: PetscScalar v,pfive = .5;
39: PetscBool flg;
40: MatFDColoring matfdcoloring = 0;
41: PetscBool fd_jacobian_coloring;
43: PetscInitialize(&argc,&argv,(char *)0,help);
44: MPI_Comm_size(PETSC_COMM_WORLD,&size);
45: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
46: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
47: h = 1.0/(n-1);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create nonlinear solver context
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: SNESCreate(PETSC_COMM_WORLD,&snes);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Create vector data structures; set function evaluation routine
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: VecCreate(PETSC_COMM_SELF,&x);
58: VecSetSizes(x,PETSC_DECIDE,n);
59: VecSetFromOptions(x);
60: VecDuplicate(x,&r);
61: VecDuplicate(x,&F);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Initialize application:
65: Store right-hand-side of PDE and exact solution
66: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: xp = 0.0;
68: for (i=0; i<n; i++) {
69: v = 6.0*xp + pow(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
70: VecSetValues(F,1,&i,&v,INSERT_VALUES);
71: xp += h;
72: }
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Evaluate initial guess; then solve nonlinear system
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: VecSet(x,pfive);
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Set Function evaluation routine
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: SNESSetFunction(snes,r,FormFunction,(void*)F);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create matrix data structures; set Jacobian evaluation routine
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&J);
88: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL,&JPrec);
89:
90: flg = PETSC_FALSE;
91: PetscOptionsGetBool(PETSC_NULL,"-snes_jacobian_default",&flg,PETSC_NULL);
92: if (flg){
93: /* Jacobian using finite differences. Slow and expensive, not take advantage of sparsity */
94: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);
95: } else {
96: /* User provided Jacobian and preconditioner(diagonal part of Jacobian) */
97: SNESSetJacobian(snes,J,JPrec,MyJacobian,0);
98: }
100: fd_jacobian_coloring = PETSC_FALSE;
101: PetscOptionsGetBool(PETSC_NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,PETSC_NULL);
102: if (fd_jacobian_coloring){
103: /* Jacobian using finite differences with matfdcoloring based on the sparse structure.
104: In this case, only three calls to FormFunction() for each Jacobian evaluation - very fast! */
105: ISColoring iscoloring;
106:
107: /* Get the data structure of J */
108: flg = PETSC_FALSE;
109: PetscOptionsGetBool(PETSC_NULL,"-my_jacobian_struct",&flg,PETSC_NULL);
110: if (flg){
111: /* use user-provided jacobian data structure */
112: MyApproxJacobianStructure(&J,PETSC_NULL);
113: } else {
114: /* use previously defined jacobian: SNESDefaultComputeJacobian() or MyJacobian() */
115: MatStructure flag;
116: SNESComputeJacobian(snes,x,&J,&J,&flag);
117: }
119: /* Create coloring context */
120: MatGetColoring(J,MATCOLORINGSL,&iscoloring);
121: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
122: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,(void*)F);
123: MatFDColoringSetFromOptions(matfdcoloring);
124: /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
125:
126: /* Use SNESDefaultComputeJacobianColor() for Jacobian evaluation */
127: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
128: ISColoringDestroy(&iscoloring);
129: }
131: SNESSetFromOptions(snes);
133: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: Solve nonlinear system
135: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: SNESSolve(snes,PETSC_NULL,x);
137: SNESGetIterationNumber(snes,&it);
138: PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Free work space. All PETSc objects should be destroyed when they
142: are no longer needed.
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: VecDestroy(&x); VecDestroy(&r);
145: VecDestroy(&F); MatDestroy(&J);
146: MatDestroy(&JPrec); SNESDestroy(&snes);
147: if (fd_jacobian_coloring){
148: MatFDColoringDestroy(&matfdcoloring);
149: }
150: PetscFinalize();
151: return 0;
152: }
153: /* ------------------------------------------------------------------- */
154: /*
155: FormInitialGuess - Forms initial approximation.
157: Input Parameters:
158: user - user-defined application context
159: X - vector
161: Output Parameter:
162: X - vector
163: */
164: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
165: {
166: PetscScalar *xx,*ff,*FF,d;
168: PetscInt i,n;
170: VecGetArray(x,&xx);
171: VecGetArray(f,&ff);
172: VecGetArray((Vec)dummy,&FF);
173: VecGetSize(x,&n);
174: d = (PetscReal)(n - 1); d = d*d;
175: ff[0] = xx[0];
176: for (i=1; i<n-1; i++) {
177: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
178: }
179: ff[n-1] = xx[n-1] - 1.0;
180: VecRestoreArray(x,&xx);
181: VecRestoreArray(f,&ff);
182: VecRestoreArray((Vec)dummy,&FF);
183: return 0;
184: }
185: /* ------------------------------------------------------------------- */
186: /*
187: MyJacobian - This routine demonstrates the use of different
188: matrices for the Jacobian and preconditioner
190: Input Parameters:
191: . snes - the SNES context
192: . x - input vector
193: . ptr - optional user-defined context, as set by SNESSetJacobian()
195: Output Parameters:
196: . A - Jacobian matrix
197: . B - different preconditioning matrix
198: . flag - flag indicating matrix structure
199: */
200: PetscErrorCode MyJacobian(SNES snes,Vec x,Mat *jac,Mat *prejac,MatStructure *flag,void *dummy)
201: {
202: PetscScalar *xx,A[3],d;
203: PetscInt i,n,j[3];
206: VecGetArray(x,&xx);
207: VecGetSize(x,&n);
208: d = (PetscReal)(n - 1); d = d*d;
210: /* Form Jacobian. Also form a different preconditioning matrix that
211: has only the diagonal elements. */
212: i = 0; A[0] = 1.0;
213: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
214: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
215: for (i=1; i<n-1; i++) {
216: j[0] = i - 1; j[1] = i; j[2] = i + 1;
217: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
218: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
219: MatSetValues(*prejac,1,&i,1,&i,&A[1],INSERT_VALUES);
220: }
221: i = n-1; A[0] = 1.0;
222: MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
223: MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
225: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
226: MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
227: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);
230: VecRestoreArray(x,&xx);
231: *flag = SAME_NONZERO_PATTERN;
232: return 0;
233: }
235: /* ------------------------------------------------------------------- */
236: /*
237: Create an approximate data structure for Jacobian matrix to be used with matcoloring
239: Input Parameters:
240: . A - dummy jacobian matrix
242: Output Parameters:
243: A - jacobian matrix with assigned non-zero structure
244: */
245: PetscErrorCode MyApproxJacobianStructure(Mat *jac,void *dummy)
246: {
247: PetscScalar zeros[3];
248: PetscInt i,n,j[3];
251: MatGetSize(*jac,&n,&n);
253: zeros[0] = zeros[1] = zeros[2] = 0.0;
254: i = 0;
255: MatSetValues(*jac,1,&i,1,&i,zeros,INSERT_VALUES);
256:
257: for (i=1; i<n-1; i++) {
258: j[0] = i - 1; j[1] = i; j[2] = i + 1;
259: MatSetValues(*jac,1,&i,3,j,zeros,INSERT_VALUES);
260: }
261: i = n-1;
262: MatSetValues(*jac,1,&i,1,&i,zeros,INSERT_VALUES);
264: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
266: return 0;
267: }