Actual source code: ex6.c

  1: /*$Id: ex6.c,v 1.71 2001/08/07 03:04:16 balay Exp $*/

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

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

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

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

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

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

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

 60:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

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

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

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

104:   SNESSetFromOptions(snes);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

222:   VecGetArray(x,&xx);
223:   VecGetSize(x,&n);
224:   d = (PetscReal)(n - 1); d = d*d;

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

241:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
242:   MatAssemblyBegin(*prejac,MAT_FINAL_ASSEMBLY);
243:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
244:   MatAssemblyEnd(*prejac,MAT_FINAL_ASSEMBLY);

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

256:    Input Parameters:
257: .  ctx - optional user-defined context, as set by PCShellSetApply()
258: .  x - input vector

260:    Output Parameter:
261: .  y - preconditioned vector
262: */
263: int MatrixFreePreconditioner(void *ctx,Vec x,Vec y)
264: {
266:   VecCopy(x,y);
267:   return 0;
268: }