Actual source code: ex7.c

petsc-3.14.1 2020-11-03
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 FormJacobianNoMatrix(SNES,Vec,Mat,Mat,void*);
  9: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 10: extern PetscErrorCode FormFunctioni(void *,PetscInt,Vec,PetscScalar *);
 11: extern PetscErrorCode OtherFunctionForDifferencing(void*,Vec,Vec);
 12: extern PetscErrorCode FormInitialGuess(SNES,Vec);
 13: extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);

 15: typedef struct {
 16:   PetscViewer viewer;
 17: } MonitorCtx;

 19: typedef struct {
 20:   PetscBool variant;
 21: } AppCtx;

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

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

 40:   /* Set up data structures */
 41:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
 42:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
 43:   VecDuplicate(x,&r);
 44:   VecDuplicate(x,&F);
 45:   VecDuplicate(x,&U);
 46:   PetscObjectSetName((PetscObject)U,"Exact Solution");

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

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

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

 64:   /* Set various routines and options */
 65:   SNESSetFunction(snes,r,FormFunction,F);
 66:   if (user.variant) {
 67:     /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
 68:     MatCreateMFFD(PETSC_COMM_WORLD,n,n,n,n,&J);
 69:     MatMFFDSetFunction(J,(PetscErrorCode (*)(void*, Vec, Vec))SNESComputeFunction,snes);
 70:     MatMFFDSetFunctioni(J,FormFunctioni);
 71:     /* Use the matrix free operator for both the Jacobian used to define the linear system and used to define the preconditioner */
 72:     /* This tests MatGetDiagonal() for MATMFFD */
 73:     PetscOptionsHasName(NULL,NULL,"-puremf",&puremf);
 74:   } else {
 75:     /* create matrix free matrix for Jacobian */
 76:     MatCreateSNESMF(snes,&J);
 77:     /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
 78:     /* note we use the same context for this function as FormFunction, the F vector */
 79:     MatMFFDSetFunction(J,OtherFunctionForDifferencing,F);
 80:   }

 82:   /* Set various routines and options */
 83:   SNESSetJacobian(snes,J,puremf ? J : B,puremf ? FormJacobianNoMatrix : FormJacobian,&user);
 84:   SNESSetFromOptions(snes);

 86:   /* Solve nonlinear system */
 87:   FormInitialGuess(snes,x);
 88:   SNESSolve(snes,NULL,x);
 89:   SNESGetIterationNumber(snes,&its);
 90:   PetscPrintf(PETSC_COMM_SELF,"number of SNES iterations = %D\n\n",its);

 92:   /* Free data structures */
 93:   VecDestroy(&x);  VecDestroy(&r);
 94:   VecDestroy(&U);  VecDestroy(&F);
 95:   MatDestroy(&J);  MatDestroy(&B);
 96:   SNESDestroy(&snes);
 97:   PetscFinalize();
 98:   return ierr;
 99: }
100: /* --------------------  Evaluate Function F(x) --------------------- */

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

109:   VecGetArrayRead(x,&xx);
110:   VecGetArray(f,&ff);
111:   VecGetArrayRead((Vec) dummy,&FF);
112:   VecGetSize(x,&n);
113:   d     = (PetscReal)(n - 1); d = d*d;
114:   ff[0] = xx[0];
115:   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];
116:   ff[n-1] = xx[n-1] - 1.0;
117:   VecRestoreArrayRead(x,&xx);
118:   VecRestoreArray(f,&ff);
119:   VecRestoreArrayRead((Vec)dummy,&FF);
120:   return 0;
121: }

123: PetscErrorCode  FormFunctioni(void *dummy,PetscInt i,Vec x,PetscScalar *s)
124: {
125:   const PetscScalar *xx,*FF;
126:   PetscScalar       d;
127:   PetscInt          n;
128:   PetscErrorCode    ierr;
129:   SNES              snes = (SNES) dummy;
130:   Vec               F;

132:   SNESGetFunction(snes,NULL,NULL,(void**)&F);
133:   VecGetArrayRead(x,&xx);
134:   VecGetArrayRead(F,&FF);
135:   VecGetSize(x,&n);
136:   d     = (PetscReal)(n - 1); d = d*d;
137:   if (i == 0) {
138:     *s = xx[0];
139:   } else if (i == n-1) {
140:     *s = xx[n-1] - 1.0;
141:   } else {
142:     *s = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
143:   }
144:   VecRestoreArrayRead(x,&xx);
145:   VecRestoreArrayRead(F,&FF);
146:   return 0;
147: }

149: /*

151:    Example function that when differenced produces the same matrix free Jacobian as FormFunction()
152:    this is provided to show how a user can provide a different function
153: */
154: PetscErrorCode  OtherFunctionForDifferencing(void *dummy,Vec x,Vec f)
155: {

158:   FormFunction(NULL,x,f,dummy);
159:   VecShift(f,1.0);
160:   return 0;
161: }

163: /* --------------------  Form initial approximation ----------------- */

165: PetscErrorCode  FormInitialGuess(SNES snes,Vec x)
166: {
168:   PetscScalar    pfive = .50;
169:   VecSet(x,pfive);
170:   return 0;
171: }
172: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
173: /*  Evaluates a matrix that is used to precondition the matrix-free
174:     jacobian. In this case, the explict preconditioner matrix is
175:     also EXACTLY the Jacobian. In general, it would be some lower
176:     order, simplified apprioximation */

178: PetscErrorCode  FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
179: {
180:   const PetscScalar *xx;
181:   PetscScalar       A[3],d;
182:   PetscInt          i,n,j[3];
183:   PetscErrorCode    ierr;
184:   AppCtx            *user = (AppCtx*) dummy;

186:   VecGetArrayRead(x,&xx);
187:   VecGetSize(x,&n);
188:   d    = (PetscReal)(n - 1); d = d*d;

190:   i    = 0; A[0] = 1.0;
191:   MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
192:   for (i=1; i<n-1; i++) {
193:     j[0] = i - 1; j[1] = i;                   j[2] = i + 1;
194:     A[0] = d;     A[1] = -2.0*d + 2.0*xx[i];  A[2] = d;
195:     MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
196:   }
197:   i     = n-1; A[0] = 1.0;
198:   MatSetValues(B,1,&i,1,&i,&A[0],INSERT_VALUES);
199:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
200:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
201:   VecRestoreArrayRead(x,&xx);

203:   if (user->variant) {
204:     MatMFFDSetBase(jac,x,NULL);
205:   }
206:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
207:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
208:   return 0;
209: }

211: PetscErrorCode  FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
212: {
213:   PetscErrorCode    ierr;
214:   AppCtx            *user = (AppCtx*) dummy;

216:   if (user->variant) {
217:     MatMFFDSetBase(jac,x,NULL);
218:   }
219:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
220:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
221:   return 0;
222: }

224: /* --------------------  User-defined monitor ----------------------- */

226: PetscErrorCode  Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
227: {
229:   MonitorCtx     *monP = (MonitorCtx*) dummy;
230:   Vec            x;
231:   MPI_Comm       comm;

233:   PetscObjectGetComm((PetscObject)snes,&comm);
234:   PetscFPrintf(comm,stdout,"iter = %D, SNES Function norm %g \n",its,(double)fnorm);
235:   SNESGetSolution(snes,&x);
236:   VecView(x,monP->viewer);
237:   return 0;
238: }


241: /*TEST

243:    test:
244:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short

246:    test:
247:       suffix: 2
248:       args: -variant -ksp_gmres_cgs_refinement_type refine_always  -snes_monitor_short
249:       output_file: output/ex7_1.out

251:    # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning
252:    test:
253:       suffix: 3
254:       args: -variant -pc_type jacobi -snes_view -ksp_monitor

256:    # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning
257:    test:
258:       suffix: 4
259:       args: -variant -pc_type jacobi -puremf  -snes_view -ksp_monitor

261: TEST*/