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: /*@C
  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: }