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: }