1 static char help[] = "Solves u`` + u^{2} = f with Newton-like methods. Using\n\
2 matrix-free techniques with user-provided explicit matrix for computing the preconditioner.\n\n";
3
4 #include <petscsnes.h>
5
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 *);
13
14 typedef struct {
15 PetscViewer viewer;
16 } MonitorCtx;
17
18 typedef struct {
19 PetscBool variant;
20 } AppCtx;
21
main(int argc,char ** argv)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;
32
33 PetscFunctionBeginUser;
34 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
35 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
36 PetscCall(PetscOptionsHasName(NULL, NULL, "-variant", &user.variant));
37 h = 1.0 / (n - 1);
38
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"));
46
47 /* create explicit matrix preconditioner */
48 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, &B));
49
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 }
58
59 /* Create nonlinear solver */
60 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
61 PetscCall(SNESSetType(snes, type));
62
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 }
80
81 /* Set various routines and options */
82 PetscCall(SNESSetJacobian(snes, J, puremf ? J : B, puremf ? FormJacobianNoMatrix : FormJacobian, &user));
83 PetscCall(SNESSetFromOptions(snes));
84
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));
90
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) --------------------- */
103
FormFunction(SNES snes,Vec x,Vec f,void * dummy)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;
109
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 }
125
FormFunctioni(void * dummy,PetscInt i,Vec x,PetscScalar * s)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;
133
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 }
152
153 /*
154
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 */
OtherFunctionForDifferencing(void * dummy,Vec x,Vec f)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 }
165
166 /* -------------------- Form initial approximation ----------------- */
167
FormInitialGuess(SNES snes,Vec x)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 matrix used to compute the preconditioner is
177 also EXACTLY the Jacobian. In general, it would be some lower
178 order, simplified apprioximation */
179
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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;
186
187 PetscFunctionBeginUser;
188 PetscCall(VecGetArrayRead(x, &xx));
189 PetscCall(VecGetSize(x, &n));
190 d = (PetscReal)(n - 1);
191 d = d * d;
192
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));
211
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 }
217
FormJacobianNoMatrix(SNES snes,Vec x,Mat jac,Mat B,void * dummy)218 PetscErrorCode FormJacobianNoMatrix(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
219 {
220 AppCtx *user = (AppCtx *)dummy;
221
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 }
228
229 /* -------------------- User-defined monitor ----------------------- */
230
Monitor(SNES snes,PetscInt its,PetscReal fnorm,void * dummy)231 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *dummy)
232 {
233 MonitorCtx *monP = (MonitorCtx *)dummy;
234 Vec x;
235 MPI_Comm comm;
236
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 }
244
245 /*TEST
246
247 test:
248 args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
249
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
254
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
260
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
265
266 TEST*/
267