Actual source code: snesmfj.c


  2: #include <petsc/private/snesimpl.h>
  3: #include <petscdm.h>
  4: #include <../src/mat/impls/mffd/mffdimpl.h>
  5: #include <petsc/private/matimpl.h>

  7: /*@C
  8:    MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which
  9:        Jacobian matrix vector products will be computed at, i.e. J(x) * a. The x is obtained
 10:        from the `SNES` object (using `SNESGetSolution()`).

 12:    Collective

 14:    Input Parameters:
 15: +   snes - the nonlinear solver context
 16: .   x - the point at which the Jacobian vector products will be performed
 17: .   jac - the matrix-free Jacobian object of `MatType` `MATMFFD`, likely obtained with `MatCreateSNESMF()`
 18: .   B - either the same as jac or another matrix type (ignored)
 19: .   flag - not relevant for matrix-free form
 20: -   dummy - the user context (ignored)

 22:    Options Database Key:
 23: .  -snes_mf - use the matrix created with `MatSNESMFCreate()` to setup the Jacobian for each new solution in the Newton process

 25:    Level: developer

 27:    Notes:
 28:       If `MatMFFDSetBase()` is ever called on jac then this routine will NO longer get
 29:     the x from the `SNES` object and `MatMFFDSetBase()` must from that point on be used to
 30:     change the base vector x.

 32:      This can be passed into `SNESSetJacobian()` as the Jacobian evaluation function argument
 33:      when using a completely matrix-free solver,
 34:      that is the B matrix is also the same matrix operator. This is used when you select
 35:      -snes_mf but rarely used directly by users. (All this routine does is call `MatAssemblyBegin/End()` on
 36:      the `Mat` jac.)

 38: .seealso: `MatMFFDGetH()`, `MatCreateSNESMF()`, `MatMFFDSetBase()`, `MatCreateMFFD()`, `MATMFFD`,
 39:           `MatMFFDSetHHistory()`, `MatMFFDSetFunctionError()`, `MatCreateMFFD()`, `SNESSetJacobian()`
 40: @*/
 41: PetscErrorCode MatMFFDComputeJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
 42: {
 43:   PetscFunctionBegin;
 44:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
 45:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat, MatAssemblyType);
 50: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat, Vec, Vec);

 52: /*@
 53:     MatSNESMFGetSNES - returns the `SNES` associated with a matrix created with `MatCreateSNESMF()`

 55:     Not Collective

 57:     Input Parameter:
 58: .   J - the matrix

 60:     Output Parameter:
 61: .   snes - the `SNES` object

 63:     Level: advanced

 65: .seealso: `Mat`, `SNES`, `MatCreateSNESMF()`
 66: @*/
 67: PetscErrorCode MatSNESMFGetSNES(Mat J, SNES *snes)
 68: {
 69:   MatMFFD j;

 71:   PetscFunctionBegin;
 72:   PetscCall(MatShellGetContext(J, &j));
 73:   *snes = (SNES)j->ctx;
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: /*
 78:    MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
 79:     base from the SNES context

 81: */
 82: static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J, MatAssemblyType mt)
 83: {
 84:   MatMFFD j;
 85:   SNES    snes;
 86:   Vec     u, f;
 87:   DM      dm;
 88:   DMSNES  dms;

 90:   PetscFunctionBegin;
 91:   PetscCall(MatShellGetContext(J, &j));
 92:   snes = (SNES)j->ctx;
 93:   PetscCall(MatAssemblyEnd_MFFD(J, mt));

 95:   PetscCall(SNESGetSolution(snes, &u));
 96:   PetscCall(SNESGetDM(snes, &dm));
 97:   PetscCall(DMGetDMSNES(dm, &dms));
 98:   if ((j->func == (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunction) && !dms->ops->computemffunction) {
 99:     PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
100:     PetscCall(MatMFFDSetBase_MFFD(J, u, f));
101:   } else {
102:     /* f value known by SNES is not correct for other differencing function */
103:     PetscCall(MatMFFDSetBase_MFFD(J, u, NULL));
104:   }
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }

108: /*
109:    MatAssemblyEnd_SNESMF_UseBase - Calls MatAssemblyEnd_MFFD() and then sets the
110:     base from the SNES context. This version will cause the base to be used for differencing
111:     even if the func is not SNESComputeFunction. See: MatSNESMFUseBase()

113: */
114: static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J, MatAssemblyType mt)
115: {
116:   MatMFFD j;
117:   SNES    snes;
118:   Vec     u, f;

120:   PetscFunctionBegin;
121:   PetscCall(MatAssemblyEnd_MFFD(J, mt));
122:   PetscCall(MatShellGetContext(J, &j));
123:   snes = (SNES)j->ctx;
124:   PetscCall(SNESGetSolution(snes, &u));
125:   PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
126:   PetscCall(MatMFFDSetBase_MFFD(J, u, f));
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: /*
131:     This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
132:   uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
133: */
134: static PetscErrorCode MatMFFDSetBase_SNESMF(Mat J, Vec U, Vec F)
135: {
136:   PetscFunctionBegin;
137:   PetscCall(MatMFFDSetBase_MFFD(J, U, F));
138:   J->ops->assemblyend = MatAssemblyEnd_MFFD;
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: static PetscErrorCode MatSNESMFSetReuseBase_SNESMF(Mat J, PetscBool use)
143: {
144:   PetscFunctionBegin;
145:   if (use) {
146:     J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase;
147:   } else {
148:     J->ops->assemblyend = MatAssemblyEnd_SNESMF;
149:   }
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: /*@
154:     MatSNESMFSetReuseBase - Causes the base vector to be used for differencing even if the function provided to `SNESSetFunction()` is not the
155:                        same as that provided to `MatMFFDSetFunction()`.

157:     Logically Collective

159:     Input Parameters:
160: +   J - the `MATMFFD` matrix
161: -   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the `MATMFFD` is
162:           not `SNESComputeFunction()`

164:     Level: advanced

166:     Note:
167:     Care must be taken when using this routine to insure that the function provided to `MatMFFDSetFunction()`, call it F_MF() is compatible with
168:     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

170:     Developer Note:
171:     This was provided for the MOOSE team who desired to have a `SNESSetFunction()` function that could change configurations (similar to variable
172:     switching) to contacts while the function provided to `MatMFFDSetFunction()` cannot. Except for the possibility of changing the configuration
173:     both functions compute the same mathematical function so the differencing makes sense.

175: .seealso: `MATMFFD`, `MatMFFDSetFunction()`, `SNESSetFunction()`, `MatCreateSNESMF()`, `MatSNESMFGetReuseBase()`
176: @*/
177: PetscErrorCode MatSNESMFSetReuseBase(Mat J, PetscBool use)
178: {
179:   PetscFunctionBegin;
181:   PetscTryMethod(J, "MatSNESMFSetReuseBase_C", (Mat, PetscBool), (J, use));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: static PetscErrorCode MatSNESMFGetReuseBase_SNESMF(Mat J, PetscBool *use)
186: {
187:   PetscFunctionBegin;
188:   if (J->ops->assemblyend == MatAssemblyEnd_SNESMF_UseBase) *use = PETSC_TRUE;
189:   else *use = PETSC_FALSE;
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: /*@
194:     MatSNESMFGetReuseBase - Determines if the base vector is to be used for differencing even if the function provided to `SNESSetFunction()` is not the
195:                        same as that provided to `MatMFFDSetFunction()`.

197:     Logically Collective

199:     Input Parameter:
200: .   J - the `MATMFFD` matrix

202:     Output Parameter:
203: .   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the `MATMFFD` is
204:           not `SNESComputeFunction()`

206:     Level: advanced

208:     Note:
209:     See `MatSNESMFSetReuseBase()`

211: .seealso: `MatSNESMFSetReuseBase()`, `MatCreateSNESMF()`, `MatSNESMFSetReuseBase()`
212: @*/
213: PetscErrorCode MatSNESMFGetReuseBase(Mat J, PetscBool *use)
214: {
215:   PetscFunctionBegin;
217:   PetscUseMethod(J, "MatSNESMFGetReuseBase_C", (Mat, PetscBool *), (J, use));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /*@
222:    MatCreateSNESMF - Creates a matrix-free matrix context for use with
223:    a `SNES` solver.  This matrix can be used as the Jacobian argument for
224:    the routine `SNESSetJacobian()`. See `MatCreateMFFD()` for details on how
225:    the finite difference computation is done.

227:    Collective

229:    Input Parameters:
230: .  snes - the `SNES` context

232:    Output Parameter:
233: .  J - the matrix-free matrix which is of type `MATMFFD`

235:    Level: advanced

237:    Notes:
238:      You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
239:      preconditioner matrix

241:      If you wish to provide a different function to do differencing on to compute the matrix-free operator than
242:      that provided to `SNESSetFunction()` then call `MatMFFDSetFunction()` with your function after this call.

244:      The difference between this routine and `MatCreateMFFD()` is that this matrix
245:      automatically gets the current base vector from the `SNES` object and not from an
246:      explicit call to `MatMFFDSetBase()`.

248:      If `MatMFFDSetBase()` is ever called on jac then this routine will NO longer get
249:      the x from the `SNES` object and `MatMFFDSetBase()` must from that point on be used to
250:      change the base vector x.

252:      Using a different function for the differencing will not work if you are using non-linear left preconditioning.

254: .seealso: `MATMFFD, `MatDestroy()`, `MatMFFDSetFunction()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`
255:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateMFFD()`,
256:           `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`, `MatSNESMFSetReuseBase()`, `MatSNESMFGetReuseBase()`
257: @*/
258: PetscErrorCode MatCreateSNESMF(SNES snes, Mat *J)
259: {
260:   PetscInt n, N;
261:   MatMFFD  mf;

263:   PetscFunctionBegin;
264:   if (snes->vec_func) {
265:     PetscCall(VecGetLocalSize(snes->vec_func, &n));
266:     PetscCall(VecGetSize(snes->vec_func, &N));
267:   } else if (snes->dm) {
268:     Vec tmp;
269:     PetscCall(DMGetGlobalVector(snes->dm, &tmp));
270:     PetscCall(VecGetLocalSize(tmp, &n));
271:     PetscCall(VecGetSize(tmp, &N));
272:     PetscCall(DMRestoreGlobalVector(snes->dm, &tmp));
273:   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
274:   PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)snes), n, n, N, N, J));
275:   PetscCall(MatShellGetContext(*J, &mf));
276:   mf->ctx = snes;

278:   if (snes->npc && snes->npcside == PC_LEFT) {
279:     PetscCall(MatMFFDSetFunction(*J, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunctionDefaultNPC, snes));
280:   } else {
281:     DM     dm;
282:     DMSNES dms;

284:     PetscCall(SNESGetDM(snes, &dm));
285:     PetscCall(DMGetDMSNES(dm, &dms));
286:     PetscCall(MatMFFDSetFunction(*J, (PetscErrorCode(*)(void *, Vec, Vec))(dms->ops->computemffunction ? SNESComputeMFFunction : SNESComputeFunction), snes));
287:   }
288:   (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;

290:   PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatMFFDSetBase_C", MatMFFDSetBase_SNESMF));
291:   PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFSetReuseBase_C", MatSNESMFSetReuseBase_SNESMF));
292:   PetscCall(PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFGetReuseBase_C", MatSNESMFGetReuseBase_SNESMF));
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }