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;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
36: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
37: PetscCall(PetscOptionsHasName(NULL, NULL, "-variant", &user.variant));
38: h = 1.0 / (n - 1);
40: /* Set up data structures */
41: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
42: PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
43: PetscCall(VecDuplicate(x, &r));
44: PetscCall(VecDuplicate(x, &F));
45: PetscCall(VecDuplicate(x, &U));
46: PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
48: /* create explicit matrix preconditioner */
49: PetscCall(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: PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
55: v = xp * xp * xp;
56: PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
57: xp += h;
58: }
60: /* Create nonlinear solver */
61: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
62: PetscCall(SNESSetType(snes, type));
64: /* Set various routines and options */
65: PetscCall(SNESSetFunction(snes, r, FormFunction, F));
66: if (user.variant) {
67: /* this approach is not normally needed, one should use the MatCreateSNESMF() below usually */
68: PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, n, n, n, n, &J));
69: PetscCall(MatMFFDSetFunction(J, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunction, snes));
70: PetscCall(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: PetscCall(PetscOptionsHasName(NULL, NULL, "-puremf", &puremf));
74: } else {
75: /* create matrix free matrix for Jacobian */
76: PetscCall(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: PetscCall(MatMFFDSetFunction(J, OtherFunctionForDifferencing, F));
80: }
82: /* Set various routines and options */
83: PetscCall(SNESSetJacobian(snes, J, puremf ? J : B, puremf ? FormJacobianNoMatrix : FormJacobian, &user));
84: PetscCall(SNESSetFromOptions(snes));
86: /* Solve nonlinear system */
87: PetscCall(FormInitialGuess(snes, x));
88: PetscCall(SNESSolve(snes, NULL, x));
89: PetscCall(SNESGetIterationNumber(snes, &its));
90: PetscCall(PetscPrintf(PETSC_COMM_SELF, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
92: /* Free data structures */
93: PetscCall(VecDestroy(&x));
94: PetscCall(VecDestroy(&r));
95: PetscCall(VecDestroy(&U));
96: PetscCall(VecDestroy(&F));
97: PetscCall(MatDestroy(&J));
98: PetscCall(MatDestroy(&B));
99: PetscCall(SNESDestroy(&snes));
100: PetscCall(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: PetscFunctionBeginUser;
112: PetscCall(VecGetArrayRead(x, &xx));
113: PetscCall(VecGetArray(f, &ff));
114: PetscCall(VecGetArrayRead((Vec)dummy, &FF));
115: PetscCall(VecGetSize(x, &n));
116: d = (PetscReal)(n - 1);
117: d = d * d;
118: ff[0] = xx[0];
119: 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];
120: ff[n - 1] = xx[n - 1] - 1.0;
121: PetscCall(VecRestoreArrayRead(x, &xx));
122: PetscCall(VecRestoreArray(f, &ff));
123: PetscCall(VecRestoreArrayRead((Vec)dummy, &FF));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: PetscErrorCode FormFunctioni(void *dummy, PetscInt i, Vec x, PetscScalar *s)
128: {
129: const PetscScalar *xx, *FF;
130: PetscScalar d;
131: PetscInt n;
132: SNES snes = (SNES)dummy;
133: Vec F;
135: PetscFunctionBeginUser;
136: PetscCall(SNESGetFunction(snes, NULL, NULL, (void **)&F));
137: PetscCall(VecGetArrayRead(x, &xx));
138: PetscCall(VecGetArrayRead(F, &FF));
139: PetscCall(VecGetSize(x, &n));
140: d = (PetscReal)(n - 1);
141: d = d * d;
142: if (i == 0) {
143: *s = xx[0];
144: } else if (i == n - 1) {
145: *s = xx[n - 1] - 1.0;
146: } else {
147: *s = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];
148: }
149: PetscCall(VecRestoreArrayRead(x, &xx));
150: PetscCall(VecRestoreArrayRead(F, &FF));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*
156: Example function that when differenced produces the same matrix free Jacobian as FormFunction()
157: this is provided to show how a user can provide a different function
158: */
159: PetscErrorCode OtherFunctionForDifferencing(void *dummy, Vec x, Vec f)
160: {
161: PetscFunctionBeginUser;
162: PetscCall(FormFunction(NULL, x, f, dummy));
163: PetscCall(VecShift(f, 1.0));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /* -------------------- Form initial approximation ----------------- */
169: PetscErrorCode FormInitialGuess(SNES snes, Vec x)
170: {
171: PetscScalar pfive = .50;
172: PetscFunctionBeginUser;
173: PetscCall(VecSet(x, pfive));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
176: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
177: /* Evaluates a matrix that is used to precondition the matrix-free
178: jacobian. In this case, the explicit preconditioner matrix is
179: also EXACTLY the Jacobian. In general, it would be some lower
180: order, simplified apprioximation */
182: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
183: {
184: const PetscScalar *xx;
185: PetscScalar A[3], d;
186: PetscInt i, n, j[3];
187: AppCtx *user = (AppCtx *)dummy;
189: PetscFunctionBeginUser;
190: PetscCall(VecGetArrayRead(x, &xx));
191: PetscCall(VecGetSize(x, &n));
192: d = (PetscReal)(n - 1);
193: d = d * d;
195: i = 0;
196: A[0] = 1.0;
197: PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES));
198: for (i = 1; i < n - 1; i++) {
199: j[0] = i - 1;
200: j[1] = i;
201: j[2] = i + 1;
202: A[0] = d;
203: A[1] = -2.0 * d + 2.0 * xx[i];
204: A[2] = d;
205: PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
206: }
207: i = n - 1;
208: A[0] = 1.0;
209: PetscCall(MatSetValues(B, 1, &i, 1, &i, &A[0], INSERT_VALUES));
210: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
211: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
212: PetscCall(VecRestoreArrayRead(x, &xx));
214: if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL));
215: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
216: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
221: {
222: AppCtx *user = (AppCtx *)dummy;
224: PetscFunctionBeginUser;
225: if (user->variant) PetscCall(MatMFFDSetBase(jac, x, NULL));
226: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
227: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /* -------------------- User-defined monitor ----------------------- */
233: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy)
234: {
235: MonitorCtx *monP = (MonitorCtx *)dummy;
236: Vec x;
237: MPI_Comm comm;
239: PetscFunctionBeginUser;
240: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
241: PetscCall(PetscFPrintf(comm, stdout, "iter = %" PetscInt_FMT ", SNES Function norm %g \n", its, (double)fnorm));
242: PetscCall(SNESGetSolution(snes, &x));
243: PetscCall(VecView(x, monP->viewer));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*TEST
249: test:
250: args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
252: test:
253: suffix: 2
254: args: -variant -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
255: output_file: output/ex7_1.out
257: # uses AIJ matrix to define diagonal matrix for Jacobian preconditioning
258: test:
259: suffix: 3
260: args: -variant -pc_type jacobi -snes_view -ksp_monitor
262: # uses MATMFFD matrix to define diagonal matrix for Jacobian preconditioning
263: test:
264: suffix: 4
265: args: -variant -pc_type jacobi -puremf -snes_view -ksp_monitor
267: TEST*/