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