Actual source code: ex7.c
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
12: typedef struct {
13: PetscViewer viewer;
14: } MonitorCtx;
16: typedef struct {
17: Mat precond;
18: PetscTruth variant;
19: } AppCtx;
23: int main(int argc,char **argv)
24: {
25: SNES snes; /* SNES context */
26: const SNESType type = SNESLS; /* default nonlinear solution method */
27: Vec x,r,F,U; /* vectors */
28: Mat J,B; /* Jacobian matrix-free, explicit preconditioner */
29: MonitorCtx monP; /* monitoring context */
30: AppCtx user; /* user-defined work context */
31: PetscScalar h,xp = 0.0,v;
32: PetscInt its,n = 5,i;
35: PetscInitialize(&argc,&argv,(char *)0,help);
36: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
37: PetscOptionsHasName(PETSC_NULL,"-variant",&user.variant);
38: h = 1.0/(n-1);
40: /* Set up data structures */
41: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&monP.viewer);
42: VecCreateSeq(PETSC_COMM_SELF,n,&x);
43: PetscObjectSetName((PetscObject)x,"Approximate Solution");
44: VecDuplicate(x,&r);
45: VecDuplicate(x,&F);
46: VecDuplicate(x,&U);
47: PetscObjectSetName((PetscObject)U,"Exact Solution");
49: /* create explict matrix preconditioner */
50: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&B);
51: user.precond = B;
53: /* Store right-hand-side of PDE and exact solution */
54: for (i=0; i<n; i++) {
55: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
56: VecSetValues(F,1,&i,&v,INSERT_VALUES);
57: v= xp*xp*xp;
58: VecSetValues(U,1,&i,&v,INSERT_VALUES);
59: xp += h;
60: }
62: /* Create nonlinear solver */
63: SNESCreate(PETSC_COMM_WORLD,&snes);
64: SNESSetType(snes,type);
66: /* Set various routines and options */
67: SNESSetFunction(snes,r,FormFunction,F);
68: if (user.variant) {
69: MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
70: MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
71: } else {
72: /* create matrix free matrix for Jacobian */
73: MatCreateSNESMF(snes,&J);
74: }
76: /* Set various routines and options */
77: SNESSetJacobian(snes,J,B,FormJacobian,&user);
78: SNESMonitorSet(snes,Monitor,&monP,0);
79: SNESSetFromOptions(snes);
81: /* Solve nonlinear system */
82: FormInitialGuess(snes,x);
83: SNESSolve(snes,PETSC_NULL,x);
84: SNESGetIterationNumber(snes,&its);
85: PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
87: /* Free data structures */
88: VecDestroy(x); VecDestroy(r);
89: VecDestroy(U); VecDestroy(F);
90: MatDestroy(J); MatDestroy(B);
91: SNESDestroy(snes);
92: PetscViewerDestroy(monP.viewer);
93: PetscFinalize();
95: return 0;
96: }
97: /* -------------------- Evaluate Function F(x) --------------------- */
99: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
100: {
101: PetscScalar *xx,*ff,*FF,d;
102: PetscInt i,n;
105: VecGetArray(x,&xx);
106: VecGetArray(f,&ff);
107: VecGetArray((Vec) dummy,&FF);
108: VecGetSize(x,&n);
109: d = (PetscReal)(n - 1); d = d*d;
110: ff[0] = xx[0];
111: for (i=1; i<n-1; i++) {
112: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
113: }
114: ff[n-1] = xx[n-1] - 1.0;
115: VecRestoreArray(x,&xx);
116: VecRestoreArray(f,&ff);
117: VecRestoreArray((Vec)dummy,&FF);
118: return 0;
119: }
120: /* -------------------- Form initial approximation ----------------- */
124: PetscErrorCode FormInitialGuess(SNES snes,Vec x)
125: {
126: PetscErrorCode ierr;
127: PetscScalar pfive = .50;
128: VecSet(x,pfive);
129: return 0;
130: }
133: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
134: /* Evaluates a matrix that is used to precondition the matrix-free
135: jacobian. In this case, the explict preconditioner matrix is
136: also EXACTLY the Jacobian. In general, it would be some lower
137: order, simplified apprioximation */
139: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)
140: {
141: PetscScalar *xx,A[3],d;
142: PetscInt i,n,j[3],iter;
144: AppCtx *user = (AppCtx*) dummy;
146: SNESGetIterationNumber(snes,&iter);
148: if (iter%2 ==0) { /* Compute new preconditioner matrix */
149: PetscPrintf(PETSC_COMM_SELF,"iter=%D, computing new preconditioning matrix\n",iter+1);
150: *B = user->precond;
151: VecGetArray(x,&xx);
152: VecGetSize(x,&n);
153: d = (PetscReal)(n - 1); d = d*d;
155: i = 0; A[0] = 1.0;
156: MatSetValues(*B,1,&i,1,&i,&A[0],INSERT_VALUES);
157: for (i=1; i<n-1; i++) {
158: j[0] = i - 1; j[1] = i; j[2] = i + 1;
159: A[0] = d; A[1] = -2.0*d + 2.0*xx[i]; A[2] = d;
160: MatSetValues(*B,1,&i,3,j,A,INSERT_VALUES);
161: }
162: i = n-1; A[0] = 1.0;
163: MatSetValues(*B,1,&i,1,&i,&A[0],INSERT_VALUES);
164: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
165: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
166: VecRestoreArray(x,&xx);
167: *flag = SAME_NONZERO_PATTERN;
168: } else { /* reuse preconditioner from last iteration */
169: PetscPrintf(PETSC_COMM_SELF,"iter=%D, using old preconditioning matrix\n",iter+1);
170: *flag = SAME_PRECONDITIONER;
171: }
172: if (user->variant) {
173: MatMFFDSetBase(*jac,x,PETSC_NULL);
174: }
175: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
176: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
177: return 0;
178: }
179: /* -------------------- User-defined monitor ----------------------- */
183: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
184: {
186: MonitorCtx *monP = (MonitorCtx*) dummy;
187: Vec x;
188: MPI_Comm comm;
190: PetscObjectGetComm((PetscObject)snes,&comm);
191: PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %G \n",its,fnorm);
192: SNESGetSolution(snes,&x);
193: VecView(x,monP->viewer);
194: return 0;
195: }