Actual source code: ex6.c

  2: static char help[] = "u`` + u^{2} = f. Different matrices for the Jacobian and the preconditioner.\n\
  3: Demonstrates the use of matrix-free Newton-Krylov methods in conjunction\n\
  4: with a user-provided preconditioner.  Input arguments are:\n\
  5:    -snes_mf : Use matrix-free Newton methods\n\
  6:    -user_precond : Employ a user-defined preconditioner.  Used only with\n\
  7:                    matrix-free methods in this example.\n\n";

  9: /*T
 10:    Concepts: SNES^different matrices for the Jacobian and preconditioner;
 11:    Concepts: SNES^matrix-free methods
 12:    Concepts: SNES^user-provided preconditioner;
 13:    Concepts: matrix-free methods
 14:    Concepts: user-provided preconditioner;
 15:    Processors: 1
 16: T*/

 18: /* 
 19:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 20:    file automatically includes:
 21:      petsc.h       - base PETSc routines   petscvec.h - vectors
 22:      petscsys.h    - system routines       petscmat.h - matrices
 23:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 24:      petscviewer.h - viewers               petscpc.h  - preconditioners
 25:      petscksp.h   - linear solvers
 26: */
 27:  #include petscsnes.h

 29: /* 
 30:    User-defined routines
 31: */
 32: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 33: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 34: PetscErrorCode MatrixFreePreconditioner(void*,Vec,Vec);

 36: int main(int argc,char **argv)
 37: {
 38:   SNES           snes;                /* SNES context */
 39:   KSP            ksp;                /* KSP context */
 40:   PC             pc;                  /* PC context */
 41:   Vec            x,r,F;               /* vectors */
 42:   Mat            J,JPrec;             /* Jacobian,preconditioner matrices */
 44:   PetscInt       it,n = 5,i;
 45:   PetscMPIInt    size;
 46:   PetscInt       *Shistit = 0,Khistl = 200,Shistl = 10;
 47:   PetscReal      h,xp = 0.0,*Khist = 0,*Shist = 0;
 48:   PetscScalar    v,pfive = .5;
 49:   PetscTruth     flg;

 51:   PetscInitialize(&argc,&argv,(char *)0,help);
 52:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 53:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");
 54:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 55:   h = 1.0/(n-1);

 57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58:      Create nonlinear solver context
 59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 61:   SNESCreate(PETSC_COMM_WORLD,&snes);

 63:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 64:      Create vector data structures; set function evaluation routine
 65:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 67:   VecCreate(PETSC_COMM_SELF,&x);
 68:   VecSetSizes(x,PETSC_DECIDE,n);
 69:   VecSetFromOptions(x);
 70:   VecDuplicate(x,&r);
 71:   VecDuplicate(x,&F);

 73:   SNESSetFunction(snes,r,FormFunction,(void*)F);

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:      Create matrix data structures; set Jacobian evaluation routine
 77:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 79:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,PETSC_NULL,&J);
 80:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,1,PETSC_NULL,&JPrec);

 82:   /*
 83:      Note that in this case we create separate matrices for the Jacobian
 84:      and preconditioner matrix.  Both of these are computed in the
 85:      routine FormJacobian()
 86:   */
 87:   SNESSetJacobian(snes,J,JPrec,FormJacobian,0);

 89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:      Customize nonlinear solver; set runtime options
 91:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 93:   /* Set preconditioner for matrix-free method */
 94:   PetscOptionsHasName(PETSC_NULL,"-snes_mf",&flg);
 95:   if (flg) {
 96:     SNESGetKSP(snes,&ksp);
 97:     KSPGetPC(ksp,&pc);
 98:     PetscOptionsHasName(PETSC_NULL,"-user_precond",&flg);
 99:     if (flg) { /* user-defined precond */
100:       PCSetType(pc,PCSHELL);
101:       PCShellSetApply(pc,MatrixFreePreconditioner);
102:     } else {PCSetType(pc,PCNONE);}
103:   }

105:   SNESSetFromOptions(snes);

107:   /*
108:      Save all the linear residuals for all the Newton steps; this enables us
109:      to retain complete convergence history for printing after the conclusion
110:      of SNESSolve().  Alternatively, one could use the monitoring options
111:            -snes_monitor -ksp_monitor
112:      to see this information during the solver's execution; however, such
113:      output during the run distorts performance evaluation data.  So, the
114:      following is a good option when monitoring code performance, for example
115:      when using -log_summary.
116:   */
117:   PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
118:   if (flg) {
119:     SNESGetKSP(snes,&ksp);
120:     PetscMalloc(Khistl*sizeof(PetscReal),&Khist);
121:     KSPSetResidualHistory(ksp,Khist,Khistl,PETSC_FALSE);
122:     PetscMalloc(Shistl*sizeof(PetscReal),&Shist);
123:     PetscMalloc(Shistl*sizeof(PetscInt),&Shistit);
124:     SNESSetConvergenceHistory(snes,Shist,Shistit,Shistl,PETSC_FALSE);
125:   }

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Initialize application:
129:      Store right-hand-side of PDE and exact solution
130:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   xp = 0.0;
133:   for (i=0; i<n; i++) {
134:     v = 6.0*xp + pow(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
135:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
136:     xp += h;
137:   }

139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:      Evaluate initial guess; then solve nonlinear system
141:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

143:   VecSet(x,pfive);
144:   SNESSolve(snes,PETSC_NULL,x);
145:   SNESGetIterationNumber(snes,&it);
146:   PetscPrintf(PETSC_COMM_SELF,"Newton iterations = %D\n\n",it);

148:   PetscOptionsHasName(PETSC_NULL,"-rhistory",&flg);
149:   if (flg) {
150:     KSPGetResidualHistory(ksp,PETSC_NULL,&Khistl);
151:     PetscRealView(Khistl,Khist,PETSC_VIEWER_STDOUT_SELF);
152:     PetscFree(Khist);
153:     SNESGetConvergenceHistory(snes,PETSC_NULL,PETSC_NULL,&Shistl);
154:     PetscRealView(Shistl,Shist,PETSC_VIEWER_STDOUT_SELF);
155:     PetscIntView(Shistl,Shistit,PETSC_VIEWER_STDOUT_SELF);
156:     PetscFree(Shist);
157:     PetscFree(Shistit);
158:   }

160:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161:      Free work space.  All PETSc objects should be destroyed when they
162:      are no longer needed.
163:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

165:   VecDestroy(x);     VecDestroy(r);
166:   VecDestroy(F);     MatDestroy(J);
167:   MatDestroy(JPrec); SNESDestroy(snes);
168:   PetscFinalize();

170:   return 0;
171: }
172: /* ------------------------------------------------------------------- */
173: /* 
174:    FormInitialGuess - Forms initial approximation.

176:    Input Parameters:
177:    user - user-defined application context
178:    X - vector

180:    Output Parameter:
181:    X - vector
182:  */
183: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy)
184: {
185:   PetscScalar    *xx,*ff,*FF,d;
187:   PetscInt       i,n;

189:   VecGetArray(x,&xx);
190:   VecGetArray(f,&ff);
191:   VecGetArray((Vec)dummy,&FF);
192:   VecGetSize(x,&n);
193:   d = (PetscReal)(n - 1); d = d*d;
194:   ff[0]   = xx[0];
195:   for (i=1; i<n-1; i++) {
196:     ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
197:   }
198:   ff[n-1] = xx[n-1] - 1.0;
199:   VecRestoreArray(x,&xx);
200:   VecRestoreArray(f,&ff);
201:   VecRestoreArray((Vec)dummy,&FF);
202:   return 0;
203: }
204: /* ------------------------------------------------------------------- */
205: /*
206:    FormJacobian - This routine demonstrates the use of different
207:    matrices for the Jacobian and preconditioner

209:    Input Parameters:
210: .  snes - the SNES context
211: .  x - input vector
212: .  ptr - optional user-defined context, as set by SNESSetJacobian()

214:    Output Parameters:
215: .  A - Jacobian matrix
216: .  B - different preconditioning matrix
217: .  flag - flag indicating matrix structure
218: */
219: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *prejac,MatStructure *flag,void *dummy)
220: {
221:   PetscScalar    *xx,A[3],d;
222:   PetscInt       i,n,j[3];

225:   VecGetArray(x,&xx);
226:   VecGetSize(x,&n);
227:   d = (PetscReal)(n - 1); d = d*d;

229:   /* Form Jacobian.  Also form a different preconditioning matrix that 
230:      has only the diagonal elements. */
231:   i = 0; A[0] = 1.0;
232:   MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
233:   MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);
234:   for (i=1; i<n-1; i++) {
235:     j[0] = i - 1; j[1] = i;                   j[2] = i + 1;
236:     A[0] = d;     A[1] = -2.0*d + 2.0*xx[i];  A[2] = d;
237:     MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
238:     MatSetValues(*prejac,1,&i,1,&i,&A[1],INSERT_VALUES);
239:   }
240:   i = n-1; A[0] = 1.0;
241:   MatSetValues(*jac,1,&i,1,&i,&A[0],INSERT_VALUES);
242:   MatSetValues(*prejac,1,&i,1,&i,&A[0],INSERT_VALUES);

244:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
245:   MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
246:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
247:   MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);

249:   VecRestoreArray(x,&xx);
250:   *flag = SAME_NONZERO_PATTERN;
251:   return 0;
252: }
253: /* ------------------------------------------------------------------- */
254: /*
255:    MatrixFreePreconditioner - This routine demonstrates the use of a
256:    user-provided preconditioner.  This code implements just the null
257:    preconditioner, which of course is not recommended for general use.

259:    Input Parameters:
260: .  ctx - optional user-defined context, as set by PCShellSetContext()
261: .  x - input vector

263:    Output Parameter:
264: .  y - preconditioned vector
265: */
266: PetscErrorCode MatrixFreePreconditioner(void *ctx,Vec x,Vec y)
267: {
269:   VecCopy(x,y);
270:   return 0;
271: }