Actual source code: ex7.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\
3: matrix-free techniques with user-provided explicit preconditioner matrix.\n\n";
5: #include <petscsnes.h>
7: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
8: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
9: extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec);
10: extern PetscErrorCode FormInitialGuess(SNES,Vec);
11: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
13: typedef struct {
14: PetscViewer viewer;
15: } MonitorCtx;
17: typedef struct {
18: Mat precond;
19: PetscBool variant;
20: } AppCtx;
24: int main(int argc,char **argv)
25: {
26: SNES snes; /* SNES context */
27: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
28: Vec x,r,F,U; /* vectors */
29: Mat J,B; /* Jacobian matrix-free, explicit preconditioner */
30: MonitorCtx monP; /* monitoring context */
31: AppCtx user; /* user-defined work context */
32: PetscScalar h,xp = 0.0,v;
33: PetscInt its,n = 5,i;
36: PetscInitialize(&argc,&argv,(char*)0,help);
37: PetscOptionsGetInt(NULL,"-n",&n,NULL);
38: PetscOptionsHasName(NULL,"-variant",&user.variant);
39: h = 1.0/(n-1);
41: /* Set up data structures */
42: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&monP.viewer);
43: VecCreateSeq(PETSC_COMM_SELF,n,&x);
44: PetscObjectSetName((PetscObject)x,"Approximate Solution");
45: VecDuplicate(x,&r);
46: VecDuplicate(x,&F);
47: VecDuplicate(x,&U);
48: PetscObjectSetName((PetscObject)U,"Exact Solution");
50: /* create explict matrix preconditioner */
51: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B);
52: user.precond = B;
54: /* Store right-hand-side of PDE and exact solution */
55: for (i=0; i<n; i++) {
56: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
57: VecSetValues(F,1,&i,&v,INSERT_VALUES);
58: v = xp*xp*xp;
59: VecSetValues(U,1,&i,&v,INSERT_VALUES);
60: xp += h;
61: }
63: /* Create nonlinear solver */
64: SNESCreate(PETSC_COMM_WORLD,&snes);
65: SNESSetType(snes,type);
67: /* Set various routines and options */
68: SNESSetFunction(snes,r,FormFunction,F);
69: if (user.variant) {
70: /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
71: MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
72: MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
73: } else {
74: /* create matrix free matrix for Jacobian */
75: MatCreateSNESMF(snes,&J);
76: /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
77: /* note we use the same context for this function as FormFunction, the F vector */
78: MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);
79: }
81: /* Set various routines and options */
82: SNESSetJacobian(snes,J,B,FormJacobian,&user);
83: SNESMonitorSet(snes,Monitor,&monP,0);
84: SNESSetFromOptions(snes);
86: /* Solve nonlinear system */
87: FormInitialGuess(snes,x);
88: SNESSolve(snes,NULL,x);
89: SNESGetIterationNumber(snes,&its);
90: PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);
92: /* Free data structures */
93: VecDestroy(&x); VecDestroy(&r);
94: VecDestroy(&U); VecDestroy(&F);
95: MatDestroy(&J); MatDestroy(&B);
96: SNESDestroy(&snes);
97: PetscViewerDestroy(&monP.viewer);
98: PetscFinalize();
100: return 0;
101: }
102: /* -------------------- Evaluate Function F(x) --------------------- */
104: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
105: {
106: PetscScalar *xx,*ff,*FF,d;
107: PetscInt i,n;
110: VecGetArray(x,&xx);
111: VecGetArray(f,&ff);
112: VecGetArray((Vec) dummy,&FF);
113: VecGetSize(x,&n);
114: d = (PetscReal)(n - 1); d = d*d;
115: ff[0] = xx[0];
116: for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
117: ff[n-1] = xx[n-1] - 1.0;
118: VecRestoreArray(x,&xx);
119: VecRestoreArray(f,&ff);
120: VecRestoreArray((Vec)dummy,&FF);
121: return 0;
122: }
124: /*
126: Example function that when differenced produces the same matrix free Jacobian as FormFunction()
127: this is provided to show how a user can provide a different function
128: */
129: PetscErrorCode OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
130: {
133: FormFunction(NULL,x,f,dummy);
134: VecShift(f,1.0);
135: return 0;
136: }
138: /* -------------------- Form initial approximation ----------------- */
142: PetscErrorCode FormInitialGuess(SNES snes,Vec x)
143: {
145: PetscScalar pfive = .50;
146: VecSet(x,pfive);
147: return 0;
148: }
151: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
152: /* Evaluates a matrix that is used to precondition the matrix-free
153: jacobian. In this case, the explict preconditioner matrix is
154: also EXACTLY the Jacobian. In general, it would be some lower
155: order, simplified apprioximation */
157: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
158: {
159: PetscScalar *xx,A[3],d;
160: PetscInt i,n,j[3],iter;
162: AppCtx *user = (AppCtx*) dummy;
164: SNESGetIterationNumber(snes,&iter);
166: if (iter%2 ==0) { /* Compute new preconditioner matrix */
167: PetscPrintf(PETSC_COMM_SELF,"iter=%D, computing new preconditioning matrix\n",iter+1);
168: *B = user->precond;
169: VecGetArray(x,&xx);
170: VecGetSize(x,&n);
171: d = (PetscReal)(n - 1); d = d*d;
173: i = 0; A[0] = 1.0;
174: MatSetValues(*B,1,&i,1,&i,&A[0],INSERT_VALUES);
175: for (i=1; i<n-1; i++) {
176: j[0] = i - 1; j[1] = i; j[2] = i + 1;
177: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
178: MatSetValues(*B,1,&i,3,j,A,INSERT_VALUES);
179: }
180: i = n-1; A[0] = 1.0;
181: MatSetValues(*B,1,&i,1,&i,&A[0],INSERT_VALUES);
182: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
183: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
184: VecRestoreArray(x,&xx);
185: *flag = SAME_NONZERO_PATTERN;
186: } else { /* reuse preconditioner from last iteration */
187: PetscPrintf(PETSC_COMM_SELF,"iter=%D, using old preconditioning matrix\n",iter+1);
188: *flag = SAME_PRECONDITIONER;
189: }
190: if (user->variant) {
191: MatMFFDSetBase(*jac,x,NULL);
192: }
193: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
194: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
195: return 0;
196: }
197: /* -------------------- User-defined monitor ----------------------- */
201: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
202: {
204: MonitorCtx *monP = (MonitorCtx*) dummy;
205: Vec x;
206: MPI_Comm comm;
208: PetscObjectGetComm((PetscObject)snes,&comm);
209: PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %G \n",its,fnorm);
210: SNESGetSolution(snes,&x);
211: VecView(x,monP->viewer);
212: return 0;
213: }