1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petscdm.h> /*I "petscdm.h" I*/
3 #include <../src/mat/impls/mffd/mffdimpl.h>
4 #include <petsc/private/matimpl.h>
5
6 /*@
7 MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which
8 Jacobian matrix-vector products will be computed at, i.e. J(x) * a. The x is obtained
9 from the `SNES` object (using `SNESGetSolution()`).
10
11 Collective
12
13 Input Parameters:
14 + snes - the nonlinear solver context
15 . x - the point at which the Jacobian-vector products will be performed
16 . jac - the matrix-free Jacobian object of `MatType` `MATMFFD`, likely obtained with `MatCreateSNESMF()`
17 . B - either the same as `jac` or another matrix type (ignored)
18 - dummy - the user context (ignored)
19
20 Options Database Key:
21 . -snes_mf - use the matrix created with `MatSNESMFCreate()` to setup the Jacobian for each new solution in the Newton process
22
23 Level: developer
24
25 Notes:
26 If `MatMFFDSetBase()` is ever called on `jac` then this routine will NO longer get
27 the `x` from the `SNES` object and `MatMFFDSetBase()` must from that point on be used to
28 change the base vector `x`.
29
30 This can be passed into `SNESSetJacobian()` as the Jacobian evaluation function argument
31 when using a completely matrix-free solver,
32 that is the B matrix is also the same matrix operator. This is used when you select
33 -snes_mf but rarely used directly by users. (All this routine does is call `MatAssemblyBegin/End()` on
34 the `Mat` `jac`.)
35
36 .seealso: [](ch_snes), `MatMFFDGetH()`, `MatCreateSNESMF()`, `MatMFFDSetBase()`, `MatCreateMFFD()`, `MATMFFD`,
37 `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`, `SNESSetJacobian()`
38 @*/
MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)39 PetscErrorCode MatMFFDComputeJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
40 {
41 PetscFunctionBegin;
42 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
43 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
47 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode MatAssemblyEnd_MFFD(Mat, MatAssemblyType);
48 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode MatMFFDSetBase_MFFD(Mat, Vec, Vec);
49
50 /*@
51 MatSNESMFGetSNES - returns the `SNES` associated with a matrix created with `MatCreateSNESMF()`
52
53 Not Collective
54
55 Input Parameter:
56 . J - the matrix
57
58 Output Parameter:
59 . snes - the `SNES` object
60
61 Level: advanced
62
63 .seealso: [](ch_snes), `Mat`, `SNES`, `MatCreateSNESMF()`
64 @*/
MatSNESMFGetSNES(Mat J,SNES * snes)65 PetscErrorCode MatSNESMFGetSNES(Mat J, SNES *snes)
66 {
67 MatMFFD j;
68
69 PetscFunctionBegin;
70 PetscCall(MatShellGetContext(J, &j));
71 *snes = (SNES)j->ctx;
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
75 /*
76 MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
77 base from the SNES context
78
79 */
MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)80 static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J, MatAssemblyType mt)
81 {
82 MatMFFD j;
83 SNES snes;
84 Vec u, f;
85 DM dm;
86 DMSNES dms;
87
88 PetscFunctionBegin;
89 PetscCall(MatShellGetContext(J, &j));
90 snes = (SNES)j->ctx;
91 PetscCall(MatAssemblyEnd_MFFD(J, mt));
92
93 PetscCall(SNESGetSolution(snes, &u));
94 PetscCall(SNESGetDM(snes, &dm));
95 PetscCall(DMGetDMSNES(dm, &dms));
96 if ((j->func == (PetscErrorCode (*)(void *, Vec, Vec))SNESComputeFunction) && !dms->ops->computemffunction) {
97 PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
98 PetscCall(MatMFFDSetBase_MFFD(J, u, f));
99 } else {
100 /* f value known by SNES is not correct for other differencing function */
101 PetscCall(MatMFFDSetBase_MFFD(J, u, NULL));
102 }
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
106 /*
107 MatAssemblyEnd_SNESMF_UseBase - Calls MatAssemblyEnd_MFFD() and then sets the
108 base from the SNES context. This version will cause the base to be used for differencing
109 even if the func is not SNESComputeFunction. See: MatSNESMFUseBase()
110
111 */
MatAssemblyEnd_SNESMF_UseBase(Mat J,MatAssemblyType mt)112 static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J, MatAssemblyType mt)
113 {
114 MatMFFD j;
115 SNES snes;
116 Vec u, f;
117
118 PetscFunctionBegin;
119 PetscCall(MatAssemblyEnd_MFFD(J, mt));
120 PetscCall(MatShellGetContext(J, &j));
121 snes = (SNES)j->ctx;
122 PetscCall(SNESGetSolution(snes, &u));
123 PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
124 PetscCall(MatMFFDSetBase_MFFD(J, u, f));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
128 /*
129 This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
130 uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
131 */
MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)132 static PetscErrorCode MatMFFDSetBase_SNESMF(Mat J, Vec U, Vec F)
133 {
134 PetscFunctionBegin;
135 PetscCall(MatMFFDSetBase_MFFD(J, U, F));
136 J->ops->assemblyend = MatAssemblyEnd_MFFD;
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
MatSNESMFSetReuseBase_SNESMF(Mat J,PetscBool use)140 static PetscErrorCode MatSNESMFSetReuseBase_SNESMF(Mat J, PetscBool use)
141 {
142 PetscFunctionBegin;
143 if (use) {
144 J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase;
145 } else {
146 J->ops->assemblyend = MatAssemblyEnd_SNESMF;
147 }
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
151 /*@
152 MatSNESMFSetReuseBase - Causes the base vector to be used for differencing even if the function provided to `SNESSetFunction()` is not the
153 same as that provided to `MatMFFDSetFunction()`.
154
155 Logically Collective
156
157 Input Parameters:
158 + J - the `MATMFFD` matrix
159 - use - if true always reuse the base vector instead of recomputing f(u) even if the function in the `MATMFFD` is
160 not `SNESComputeFunction()`
161
162 Level: advanced
163
164 Note:
165 Care must be taken when using this routine to insure that the function provided to `MatMFFDSetFunction()`, call it F_MF() is compatible with
166 with that provided to `SNESSetFunction()`, call it F_SNES(). That is, (F_MF(u + h*d) - F_SNES(u))/h has to approximate the derivative
167
168 Developer Notes:
169 This was provided for the MOOSE team who desired to have a `SNESSetFunction()` function that could change configurations (similar to variable
170 switching) to contacts while the function provided to `MatMFFDSetFunction()` cannot. Except for the possibility of changing the configuration
171 both functions compute the same mathematical function so the differencing makes sense.
172
173 .seealso: [](ch_snes), `SNES`, `MATMFFD`, `MatMFFDSetFunction()`, `SNESSetFunction()`, `MatCreateSNESMF()`, `MatSNESMFGetReuseBase()`
174 @*/
MatSNESMFSetReuseBase(Mat J,PetscBool use)175 PetscErrorCode MatSNESMFSetReuseBase(Mat J, PetscBool use)
176 {
177 PetscFunctionBegin;
178 PetscValidHeaderSpecific(J, MAT_CLASSID, 1);
179 PetscTryMethod(J, "MatSNESMFSetReuseBase_C", (Mat, PetscBool), (J, use));
180 PetscFunctionReturn(PETSC_SUCCESS);
181 }
182
MatSNESMFGetReuseBase_SNESMF(Mat J,PetscBool * use)183 static PetscErrorCode MatSNESMFGetReuseBase_SNESMF(Mat J, PetscBool *use)
184 {
185 PetscFunctionBegin;
186 if (J->ops->assemblyend == MatAssemblyEnd_SNESMF_UseBase) *use = PETSC_TRUE;
187 else *use = PETSC_FALSE;
188 PetscFunctionReturn(PETSC_SUCCESS);
189 }
190
191 /*@
192 MatSNESMFGetReuseBase - Determines if the base vector is to be used for differencing even if the function provided to `SNESSetFunction()` is not the
193 same as that provided to `MatMFFDSetFunction()`.
194
195 Logically Collective
196
197 Input Parameter:
198 . J - the `MATMFFD` matrix
199
200 Output Parameter:
201 . use - if true always reuse the base vector instead of recomputing f(u) even if the function in the `MATMFFD` is
202 not `SNESComputeFunction()`
203
204 Level: advanced
205
206 Note:
207 See `MatSNESMFSetReuseBase()`
208
209 .seealso: [](ch_snes), `Mat`, `SNES`, `MatSNESMFSetReuseBase()`, `MatCreateSNESMF()`
210 @*/
MatSNESMFGetReuseBase(Mat J,PetscBool * use)211 PetscErrorCode MatSNESMFGetReuseBase(Mat J, PetscBool *use)
212 {
213 PetscFunctionBegin;
214 PetscValidHeaderSpecific(J, MAT_CLASSID, 1);
215 PetscUseMethod(J, "MatSNESMFGetReuseBase_C", (Mat, PetscBool *), (J, use));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
219 /*@
220 MatCreateSNESMF - Creates a finite differencing based matrix-free matrix context for use with
221 a `SNES` solver. This matrix can be used as the Jacobian argument for
222 the routine `SNESSetJacobian()`. See `MatCreateMFFD()` for details on how
223 the finite difference computation is done.
224
225 Collective
226
227 Input Parameters:
228 . snes - the `SNES` context
229
230 Output Parameter:
231 . J - the matrix-free matrix which is of type `MATMFFD`
232
233 Level: advanced
234
235 Notes:
236 You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are not using a different
237 matrix to construct the preconditioner.
238
239 If you wish to provide a different function to do differencing on to compute the matrix-free operator than
240 that provided to `SNESSetFunction()` then call `MatMFFDSetFunction()` with your function after this call.
241
242 The difference between this routine and `MatCreateMFFD()` is that this matrix
243 automatically gets the current base vector from the `SNES` object and not from an
244 explicit call to `MatMFFDSetBase()`.
245
246 If `MatMFFDSetBase()` is ever called on `jac` then this routine will NO longer get
247 the x from the `SNES` object and `MatMFFDSetBase()` must from that point on be used to
248 change the base vector `x`.
249
250 Using a different function for the differencing will not work if you are using non-linear left preconditioning.
251
252 This uses finite-differencing to apply the operator. To create a matrix-free `Mat` whose matrix-vector operator you
253 provide with your own function use `MatCreateShell()`.
254
255 Developer Note:
256 This function should really be called `MatCreateSNESMFFD()` in correspondence to `MatCreateMFFD()` to clearly indicate
257 that this is for using finite differences to apply the operator matrix-free.
258
259 .seealso: [](ch_snes), `SNES`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunction()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`
260 `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateMFFD()`, `MatCreateShell()`,
261 `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`, `MatSNESMFSetReuseBase()`, `MatSNESMFGetReuseBase()`
262 @*/
MatCreateSNESMF(SNES snes,Mat * J)263 PetscErrorCode MatCreateSNESMF(SNES snes, Mat *J)
264 {
265 PetscInt n, N;
266 MatMFFD mf;
267
268 PetscFunctionBegin;
269 if (snes->vec_func) {
270 PetscCall(VecGetLocalSize(snes->vec_func, &n));
271 PetscCall(VecGetSize(snes->vec_func, &N));
272 } else if (snes->dm) {
273 Vec tmp;
274 PetscCall(DMGetGlobalVector(snes->dm, &tmp));
275 PetscCall(VecGetLocalSize(tmp, &n));
276 PetscCall(VecGetSize(tmp, &N));
277 PetscCall(DMRestoreGlobalVector(snes->dm, &tmp));
278 } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
279 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)snes), n, n, N, N, J));
280 PetscCall(MatShellGetContext(*J, &mf));
281 mf->ctx = snes;
282
283 if (snes->npc && snes->npcside == PC_LEFT) {
284 PetscCall(MatMFFDSetFunction(*J, (PetscErrorCode (*)(void *, Vec, Vec))SNESComputeFunctionDefaultNPC, snes));
285 } else {
286 DM dm;
287 DMSNES dms;
288
289 PetscCall(SNESGetDM(snes, &dm));
290 PetscCall(DMGetDMSNES(dm, &dms));
291 PetscCall(MatMFFDSetFunction(*J, (PetscErrorCode (*)(void *, Vec, Vec))(dms->ops->computemffunction ? SNESComputeMFFunction : SNESComputeFunction), snes));
292 }
293 (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;
294
295 PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatMFFDSetBase_C", MatMFFDSetBase_SNESMF));
296 PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFSetReuseBase_C", MatSNESMFSetReuseBase_SNESMF));
297 PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFGetReuseBase_C", MatSNESMFGetReuseBase_SNESMF));
298 PetscFunctionReturn(PETSC_SUCCESS);
299 }
300