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>

  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;
 32:   PetscBool   puremf = PETSC_FALSE;

 35:   PetscInitialize(&argc, &argv, (char *)0, help);
 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 explicit 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 = %" PetscInt_FMT "\n\n", its);

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

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

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

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

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

152: /*

154:    Example function that when differenced produces the same matrix free Jacobian as FormFunction()
155:    this is provided to show how a user can provide a different function
156: */
157: PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f)
158: {
159:   FormFunction(NULL, x, f, dummy);
160:   VecShift(f, 1.0);
161:   return 0;
162: }

164: /* --------------------  Form initial approximation ----------------- */

166: PetscErrorCode FormInitialGuess(SNES snes, Vec x)
167: {
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 explicit 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:   AppCtx            *user = (AppCtx *)dummy;

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

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

209:   if (user->variant) MatMFFDSetBase(jac, x, NULL);
210:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
211:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
212:   return 0;
213: }

215: PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
216: {
217:   AppCtx *user = (AppCtx *)dummy;

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

225: /* --------------------  User-defined monitor ----------------------- */

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

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

240: /*TEST

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

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

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

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

260: TEST*/