Actual source code: snesmfj.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
3: #include <../src/mat/impls/mffd/mffdimpl.h>
4: #include <petsc/private/matimpl.h>
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()`).
11: Collective
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)
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
23: Level: developer
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`.
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`.)
36: .seealso: [](ch_snes), `MatMFFDGetH()`, `MatCreateSNESMF()`, `MatMFFDSetBase()`, `MatCreateMFFD()`, `MATMFFD`,
37: `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`, `SNESSetJacobian()`
38: @*/
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: }
47: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat, MatAssemblyType);
48: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat, Vec, Vec);
50: /*@
51: MatSNESMFGetSNES - returns the `SNES` associated with a matrix created with `MatCreateSNESMF()`
53: Not Collective
55: Input Parameter:
56: . J - the matrix
58: Output Parameter:
59: . snes - the `SNES` object
61: Level: advanced
63: .seealso: [](ch_snes), `Mat`, `SNES`, `MatCreateSNESMF()`
64: @*/
65: PetscErrorCode MatSNESMFGetSNES(Mat J, SNES *snes)
66: {
67: MatMFFD j;
69: PetscFunctionBegin;
70: PetscCall(MatShellGetContext(J, &j));
71: *snes = (SNES)j->ctx;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: /*
76: MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
77: base from the SNES context
79: */
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;
88: PetscFunctionBegin;
89: PetscCall(MatShellGetContext(J, &j));
90: snes = (SNES)j->ctx;
91: PetscCall(MatAssemblyEnd_MFFD(J, mt));
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: }
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()
111: */
112: static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J, MatAssemblyType mt)
113: {
114: MatMFFD j;
115: SNES snes;
116: Vec u, f;
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: }
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: */
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: }
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: }
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()`.
155: Logically Collective
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()`
162: Level: advanced
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
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.
173: .seealso: [](ch_snes), `SNES`, `MATMFFD`, `MatMFFDSetFunction()`, `SNESSetFunction()`, `MatCreateSNESMF()`, `MatSNESMFGetReuseBase()`
174: @*/
175: PetscErrorCode MatSNESMFSetReuseBase(Mat J, PetscBool use)
176: {
177: PetscFunctionBegin;
179: PetscTryMethod(J, "MatSNESMFSetReuseBase_C", (Mat, PetscBool), (J, use));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
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: }
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()`.
195: Logically Collective
197: Input Parameter:
198: . J - the `MATMFFD` matrix
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()`
204: Level: advanced
206: Note:
207: See `MatSNESMFSetReuseBase()`
209: .seealso: [](ch_snes), `Mat`, `SNES`, `MatSNESMFSetReuseBase()`, `MatCreateSNESMF()`
210: @*/
211: PetscErrorCode MatSNESMFGetReuseBase(Mat J, PetscBool *use)
212: {
213: PetscFunctionBegin;
215: PetscUseMethod(J, "MatSNESMFGetReuseBase_C", (Mat, PetscBool *), (J, use));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
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.
225: Collective
227: Input Parameters:
228: . snes - the `SNES` context
230: Output Parameter:
231: . J - the matrix-free matrix which is of type `MATMFFD`
233: Level: advanced
235: Notes:
236: You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are not using a different
237: preconditioner matrix
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.
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()`.
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`.
250: Using a different function for the differencing will not work if you are using non-linear left preconditioning.
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()`.
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.
259: .seealso: [](ch_snes), `SNES`, `MATMFFD`, `MatDestroy()`, `MatMFFDSetFunction()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`
260: `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateMFFD()`, `MatCreateShell()`,
261: `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`, `MatSNESMFSetReuseBase()`, `MatSNESMFGetReuseBase()`
262: @*/
263: PetscErrorCode MatCreateSNESMF(SNES snes, Mat *J)
264: {
265: PetscInt n, N;
266: MatMFFD mf;
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;
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;
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;
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: }