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:    Option 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:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
 44:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
 45:   return 0;
 46: }

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

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

 54:     Not collective

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

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

 62:     Level: advanced

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

 70:   MatShellGetContext(J, &j);
 71:   *snes = (SNES)j->ctx;
 72:   return 0;
 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:   MatShellGetContext(J, &j);
 89:   snes = (SNES)j->ctx;
 90:   MatAssemblyEnd_MFFD(J, mt);

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

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

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

117:   MatAssemblyEnd_MFFD(J, mt);
118:   MatShellGetContext(J, &j);
119:   snes = (SNES)j->ctx;
120:   SNESGetSolution(snes, &u);
121:   SNESGetFunction(snes, &f, NULL, NULL);
122:   MatMFFDSetBase_MFFD(J, u, f);
123:   return 0;
124: }

126: /*
127:     This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
128:   uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
129: */
130: static PetscErrorCode MatMFFDSetBase_SNESMF(Mat J, Vec U, Vec F)
131: {
132:   MatMFFDSetBase_MFFD(J, U, F);
133:   J->ops->assemblyend = MatAssemblyEnd_MFFD;
134:   return 0;
135: }

137: static PetscErrorCode MatSNESMFSetReuseBase_SNESMF(Mat J, PetscBool use)
138: {
139:   if (use) {
140:     J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase;
141:   } else {
142:     J->ops->assemblyend = MatAssemblyEnd_SNESMF;
143:   }
144:   return 0;
145: }

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

151:     Logically Collective

153:     Input Parameters:
154: +   J - the `MATMFFD` matrix
155: -   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the `MATMFFD` is
156:           not `SNESComputeFunction()`

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

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

167:     Level: advanced

169: .seealso: `MATMFFD`, `MatMFFDSetFunction()`, `SNESSetFunction()`, `MatCreateSNESMF()`, `MatSNESMFGetReuseBase()`
170: @*/
171: PetscErrorCode MatSNESMFSetReuseBase(Mat J, PetscBool use)
172: {
174:   PetscTryMethod(J, "MatSNESMFSetReuseBase_C", (Mat, PetscBool), (J, use));
175:   return 0;
176: }

178: static PetscErrorCode MatSNESMFGetReuseBase_SNESMF(Mat J, PetscBool *use)
179: {
180:   if (J->ops->assemblyend == MatAssemblyEnd_SNESMF_UseBase) *use = PETSC_TRUE;
181:   else *use = PETSC_FALSE;
182:   return 0;
183: }

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

189:     Logically Collective

191:     Input Parameter:
192: .   J - the `MATMFFD` matrix

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

198:     Level: advanced

200:     Note:
201:     See `MatSNESMFSetReuseBase()`

203: .seealso: `MatSNESMFSetReuseBase()`, `MatCreateSNESMF()`, `MatSNESMFSetReuseBase()`
204: @*/
205: PetscErrorCode MatSNESMFGetReuseBase(Mat J, PetscBool *use)
206: {
208:   PetscUseMethod(J, "MatSNESMFGetReuseBase_C", (Mat, PetscBool *), (J, use));
209:   return 0;
210: }

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

218:    Collective

220:    Input Parameters:
221: .  snes - the `SNES` context

223:    Output Parameter:
224: .  J - the matrix-free matrix which is of type `MATMFFD`

226:    Level: advanced

228:    Notes:
229:      You can call `SNESSetJacobian()` with `MatMFFDComputeJacobian()` if you are using matrix and not a different
230:      preconditioner matrix

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

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

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

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

245: .seealso: `MATMFFD, `MatDestroy()`, `MatMFFDSetFunction()`, `MatMFFDSetFunctionError()`, `MatMFFDDSSetUmin()`
246:           `MatMFFDSetHHistory()`, `MatMFFDResetHHistory()`, `MatCreateMFFD()`,
247:           `MatMFFDGetH()`, `MatMFFDRegister()`, `MatMFFDComputeJacobian()`, `MatSNESMFSetReuseBase()`, `MatSNESMFGetReuseBase()`
248: @*/
249: PetscErrorCode MatCreateSNESMF(SNES snes, Mat *J)
250: {
251:   PetscInt n, N;
252:   MatMFFD  mf;

254:   if (snes->vec_func) {
255:     VecGetLocalSize(snes->vec_func, &n);
256:     VecGetSize(snes->vec_func, &N);
257:   } else if (snes->dm) {
258:     Vec tmp;
259:     DMGetGlobalVector(snes->dm, &tmp);
260:     VecGetLocalSize(tmp, &n);
261:     VecGetSize(tmp, &N);
262:     DMRestoreGlobalVector(snes->dm, &tmp);
263:   } else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
264:   MatCreateMFFD(PetscObjectComm((PetscObject)snes), n, n, N, N, J);
265:   MatShellGetContext(*J, &mf);
266:   mf->ctx = snes;

268:   if (snes->npc && snes->npcside == PC_LEFT) {
269:     MatMFFDSetFunction(*J, (PetscErrorCode(*)(void *, Vec, Vec))SNESComputeFunctionDefaultNPC, snes);
270:   } else {
271:     DM     dm;
272:     DMSNES dms;

274:     SNESGetDM(snes, &dm);
275:     DMGetDMSNES(dm, &dms);
276:     MatMFFDSetFunction(*J, (PetscErrorCode(*)(void *, Vec, Vec))(dms->ops->computemffunction ? SNESComputeMFFunction : SNESComputeFunction), snes);
277:   }
278:   (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;

280:   PetscObjectComposeFunction((PetscObject)*J, "MatMFFDSetBase_C", MatMFFDSetBase_SNESMF);
281:   PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFSetReuseBase_C", MatSNESMFSetReuseBase_SNESMF);
282:   PetscObjectComposeFunction((PetscObject)*J, "MatSNESMFGetReuseBase_C", MatSNESMFGetReuseBase_SNESMF);
283:   return 0;
284: }