Actual source code: ex7.c

petsc-master 2017-01-17
Report Typos and Errors

  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>

  7: extern PetscErrorCode   FormJacobian(SNES,Vec,Mat,Mat,void*);
  8: extern PetscErrorCode   FormFunction(SNES,Vec,Vec,void*);
  9: extern PetscErrorCode   OtherFunctionForDifferencing(void*,Vec,Vec);
 10: extern PetscErrorCode   FormInitialGuess(SNES,Vec);
 11: extern PetscErrorCode   Monitor(SNES,PetscInt,PetscReal,void*);

 13: typedef struct {
 14:   PetscViewer viewer;
 15: } MonitorCtx;

 17: typedef struct {
 18:   PetscBool variant;
 19: } AppCtx;

 21: int main(int argc,char **argv)
 22: {
 23:   SNES           snes;                 /* SNES context */
 24:   SNESType       type = SNESNEWTONLS;        /* default nonlinear solution method */
 25:   Vec            x,r,F,U;              /* vectors */
 26:   Mat            J,B;                  /* Jacobian matrix-free, explicit preconditioner */
 27:   MonitorCtx     monP;                 /* monitoring context */
 28:   AppCtx         user;                 /* user-defined work context */
 29:   PetscScalar    h,xp = 0.0,v;
 30:   PetscInt       its,n = 5,i;

 33:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 34:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 35:   PetscOptionsHasName(NULL,NULL,"-variant",&user.variant);
 36:   h    = 1.0/(n-1);

 38:   /* Set up data structures */
 39:   PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&monP.viewer);
 40:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
 41:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
 42:   VecDuplicate(x,&r);
 43:   VecDuplicate(x,&F);
 44:   VecDuplicate(x,&U);
 45:   PetscObjectSetName((PetscObject)U,"Exact Solution");

 47:   /* create explict matrix preconditioner */
 48:   MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&B);

 50:   /* Store right-hand-side of PDE and exact solution */
 51:   for (i=0; i<n; i++) {
 52:     v    = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
 53:     VecSetValues(F,1,&i,&v,INSERT_VALUES);
 54:     v    = xp*xp*xp;
 55:     VecSetValues(U,1,&i,&v,INSERT_VALUES);
 56:     xp  += h;
 57:   }

 59:   /* Create nonlinear solver */
 60:   SNESCreate(PETSC_COMM_WORLD,&snes);
 61:   SNESSetType(snes,type);

 63:   /* Set various routines and options */
 64:   SNESSetFunction(snes,r,FormFunction,F);
 65:   if (user.variant) {
 66:     /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
 67:     MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
 68:     MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
 69:   } else {
 70:     /* create matrix free matrix for Jacobian */
 71:     MatCreateSNESMF(snes,&J);
 72:     /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
 73:     /* note we use the same context for this function as FormFunction, the F vector */
 74:     MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);
 75:   }

 77:   /* Set various routines and options */
 78:   SNESSetJacobian(snes,J,B,FormJacobian,&user);
 79:   SNESMonitorSet(snes,Monitor,&monP,0);
 80:   SNESSetFromOptions(snes);

 82:   /* Solve nonlinear system */
 83:   FormInitialGuess(snes,x);
 84:   SNESSolve(snes,NULL,x);
 85:   SNESGetIterationNumber(snes,&its);
 86:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

 88:   /* Free data structures */
 89:   VecDestroy(&x);  VecDestroy(&r);
 90:   VecDestroy(&U);  VecDestroy(&F);
 91:   MatDestroy(&J);  MatDestroy(&B);
 92:   SNESDestroy(&snes);
 93:   PetscViewerDestroy(&monP.viewer);
 94:   PetscFinalize();
 95:   return ierr;
 96: }
 97: /* --------------------  Evaluate Function F(x) --------------------- */

 99: PetscErrorCode  FormFunction(SNES snes,Vec x,Vec f,void *dummy)
100: {
101:   const PetscScalar *xx;
102:   PetscScalar       *ff,*FF,d;
103:   PetscInt          i,n;
104:   PetscErrorCode    ierr;

106:   VecGetArrayRead(x,&xx);
107:   VecGetArray(f,&ff);
108:   VecGetArray((Vec) dummy,&FF);
109:   VecGetSize(x,&n);
110:   d     = (PetscReal)(n - 1); d = d*d;
111:   ff[0] = xx[0];
112:   for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
113:   ff[n-1] = xx[n-1] - 1.0;
114:   VecRestoreArrayRead(x,&xx);
115:   VecRestoreArray(f,&ff);
116:   VecRestoreArray((Vec)dummy,&FF);
117:   return 0;
118: }

120: /*

122:    Example function that when differenced produces the same matrix free Jacobian as FormFunction()
123:    this is provided to show how a user can provide a different function
124: */
125: PetscErrorCode  OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
126: {

129:   FormFunction(NULL,x,f,dummy);
130:   VecShift(f,1.0);
131:   return 0;
132: }

134: /* --------------------  Form initial approximation ----------------- */

136: PetscErrorCode  FormInitialGuess(SNES snes,Vec x)
137: {
139:   PetscScalar    pfive = .50;
140:   VecSet(x,pfive);
141:   return 0;
142: }
143: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
144: /*  Evaluates a matrix that is used to precondition the matrix-free
145:     jacobian. In this case, the explict preconditioner matrix is
146:     also EXACTLY the Jacobian. In general, it would be some lower
147:     order, simplified apprioximation */

149: PetscErrorCode  FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
150: {
151:   const PetscScalar *xx;
152:   PetscScalar       A[3],d;
153:   PetscInt          i,n,j[3];
154:   PetscErrorCode    ierr;
155:   AppCtx            *user = (AppCtx*) dummy;

157:   VecGetArrayRead(x,&xx);
158:   VecGetSize(x,&n);
159:   d    = (PetscReal)(n - 1); d = d*d;

161:   i    = 0; A[0] = 1.0;
162:   MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
163:   for (i=1; i<n-1; i++) {
164:     j[0] = i - 1; j[1] = i;                   j[2] = i + 1;
165:     A[0] = d;     A[1] = -2.0*d + 2.0*xx[i];  A[2] = d;
166:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
167:   }
168:   i     = n-1; A[0] = 1.0;
169:   MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
170:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
171:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
172:   VecRestoreArrayRead(x,&xx);

174:   if (user->variant) {
175:     MatMFFDSetBase(jac,x,NULL);
176:   }
177:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
178:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
179:   return 0;
180: }
181: /* --------------------  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,(double)fnorm);
192:   SNESGetSolution(snes,&x);
193:   VecView(x,monP->viewer);
194:   return 0;
195: }