Actual source code: ex7.c
1: static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\
2: matrix-free techniques with user-provided explicit preconditioner matrix.\n\n";
4: #include <petscsnes.h>
6: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
7: extern PetscErrorCode FormJacobianNoMatrix(SNES, Vec, Mat, Mat, void *);
8: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
9: extern PetscErrorCode FormFunctioni(void *, PetscInt, Vec, PetscScalar *);
10: extern PetscErrorCode OtherFunctionForDifferencing(void *, Vec, Vec);
11: extern PetscErrorCode FormInitialGuess(SNES, Vec);
12: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
14: typedef struct {
15: PetscViewer viewer;
16: } MonitorCtx;
18: typedef struct {
19: PetscBool variant;
20: } AppCtx;
22: int main(int argc, char **argv)
23: {
24: SNES snes; /* SNES context */
25: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
26: Vec x, r, F, U; /* vectors */
27: Mat J, B; /* Jacobian matrix-free, explicit preconditioner */
28: AppCtx user; /* user-defined work context */
29: PetscScalar h, xp = 0.0, v;
30: PetscInt its, n = 5, i;
31: PetscBool puremf = PETSC_FALSE;
33: PetscFunctionBeginUser;
34: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
35: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
36: PetscCall(PetscOptionsHasName(NULL, NULL, "-variant", &user.variant));
37: h = 1.0 / (n - 1);
39: /* Set up data structures */
40: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
41: PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
42: PetscCall(VecDuplicate(x, &r));
43: PetscCall(VecDuplicate(x, &F));
44: PetscCall(VecDuplicate(x, &U));
45: PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
47: /* create explicit matrix preconditioner */
48: PetscCall(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: PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
54: v = xp * xp * xp;
55: PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
56: xp += h;
57: }
59: /* Create nonlinear solver */
60: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
61: PetscCall(SNESSetType(snes, type));
63: /* Set various routines and options */
64: PetscCall(SNESSetFunction(snes, r, FormFunction, F));
65: if (user.variant) {
66: /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
67: PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, n, n, n, n, &J));
68: PetscCall(MatMFFDSetFunction(J, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunction, snes));
69: PetscCall(MatMFFDSetFunctioni(J, FormFunctioni));
70: /* Use the matrix-free operator for both the Jacobian used to define the linear system and used to define the preconditioner */
71: /* This tests MatGetDiagonal() for MATMFFD */
72: PetscCall(PetscOptionsHasName(NULL, NULL, "-puremf", &puremf));
73: } else {
74: /* create matrix-free matrix for Jacobian */
75: PetscCall(MatCreateSNESMF(snes, &J));
76: /* demonstrates differencing a different function than FormFunction() to apply a matrix operator */
77: /* note we use the same context for this function as FormFunction, the F vector */
78: PetscCall(MatMFFDSetFunction(J, OtherFunctionForDifferencing, F));
79: }
81: /* Set various routines and options */
82: PetscCall(SNESSetJacobian(snes, J, puremf ? J : B, puremf ? FormJacobianNoMatrix : FormJacobian, &user));
83: PetscCall(SNESSetFromOptions(snes));
85: /* Solve nonlinear system */
86: PetscCall(FormInitialGuess(snes, x));
87: PetscCall(SNESSolve(snes, NULL, x));
88: PetscCall(SNESGetIterationNumber(snes, &its));
89: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
91: /* Free data structures */
92: PetscCall(VecDestroy(&x));
93: PetscCall(VecDestroy(&r));
94: PetscCall(VecDestroy(&U));
95: PetscCall(VecDestroy(&F));
96: PetscCall(MatDestroy(&J));
97: PetscCall(MatDestroy(&B));
98: PetscCall(SNESDestroy(&snes));
99: PetscCall(PetscFinalize());
100: return 0;
101: }
102: /* -------------------- Evaluate Function F(x) --------------------- */
104: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *dummy)
105: {
106: const PetscScalar *xx, *FF;
107: PetscScalar *ff, d;
108: PetscInt i, n;
110: PetscFunctionBeginUser;
111: PetscCall(VecGetArrayRead(x, &xx));
112: PetscCall(VecGetArray(f, &ff));
113: PetscCall(VecGetArrayRead((Vec)dummy, &FF));
114: PetscCall(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: PetscCall(VecRestoreArrayRead(x, &xx));
121: PetscCall(VecRestoreArray(f, &ff));
122: PetscCall(VecRestoreArrayRead((Vec)dummy, &FF));
123: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBeginUser;
135: PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F));
136: PetscCall(VecGetArrayRead(x, &xx));
137: PetscCall(VecGetArrayRead(F, &FF));
138: PetscCall(VecGetSize(x, &n));
139: d = (PetscReal)(n - 1);
140: d = d * d;
141: if (i == 0) {
142: *s = xx[0];
143: } else if (i == n - 1) {
144: *s = xx[n - 1] - 1.0;
145: } else {
146: *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];
147: }
148: PetscCall(VecRestoreArrayRead(x, &xx));
149: PetscCall(VecRestoreArrayRead(F, &FF));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /*
155: Example function that when differenced produces the same matrix-free Jacobian as FormFunction()
156: this is provided to show how a user can provide a different function
157: */
158: PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f)
159: {
160: PetscFunctionBeginUser;
161: PetscCall(FormFunction(NULL, x, f, dummy));
162: PetscCall(VecShift(f, 1.0));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /* -------------------- Form initial approximation ----------------- */
168: PetscErrorCode FormInitialGuess(SNES snes, Vec x)
169: {
170: PetscFunctionBeginUser;
171: PetscCall(VecSet(x, 0.5));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
174: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
175: /* Evaluates a matrix that is used to precondition the matrix-free
176: jacobian. In this case, the explicit preconditioner matrix is
177: also EXACTLY the Jacobian. In general, it would be some lower
178: order, simplified apprioximation */
180: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
181: {
182: const PetscScalar *xx;
183: PetscScalar A[3], d;
184: PetscInt i, n, j[3];
185: AppCtx *user = (AppCtx *)dummy;
187: PetscFunctionBeginUser;
188: PetscCall(VecGetArrayRead(x, &xx));
189: PetscCall(VecGetSize(x, &n));
190: d = (PetscReal)(n - 1);
191: d = d * d;
193: i = 0;
194: A[0] = 1.0;
195: PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES));
196: for (i = 1; i < n - 1; i++) {
197: j[0] = i - 1;
198: j[1] = i;
199: j[2] = i + 1;
200: A[0] = d;
201: A[1] = -2.0 * d + 2.0 * xx[i];
202: A[2] = d;
203: PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
204: }
205: i = n - 1;
206: A[0] = 1.0;
207: PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES));
208: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
209: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
210: PetscCall(VecRestoreArrayRead(x, &xx));
212: if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL));
213: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
214: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
219: {
220: AppCtx *user = (AppCtx *)dummy;
222: PetscFunctionBeginUser;
223: if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL));
224: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
225: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: /* -------------------- User-defined monitor ----------------------- */
231: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy)
232: {
233: MonitorCtx *monP = (MonitorCtx *)dummy;
234: Vec x;
235: MPI_Comm comm;
237: PetscFunctionBeginUser;
238: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
239: PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm));
240: PetscCall(SNESGetSolution(snes, &x));
241: PetscCall(VecView(x, monP->viewer));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*TEST
247: test:
248: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
250: test:
251: suffix: 2
252: args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
253: output_file: output/ex7_1.out
255: # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning
256: test:
257: requires: !single
258: suffix: 3
259: args: -variant -pc_type jacobi -snes_view -ksp_monitor
261: # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning
262: test:
263: suffix: 4
264: args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor
266: TEST*/