Actual source code: snesmfj.c
petsc-3.6.1 2015-08-06
2: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3: #include <petscdm.h> /*I "petscdm.h" I*/
4: #include <../src/mat/impls/mffd/mffdimpl.h>
5: #include <petsc/private/matimpl.h>
9: /*@C
10: MatMFFDComputeJacobian - Tells the matrix-free Jacobian object the new location at which
11: Jacobian matrix vector products will be computed at, i.e. J(x) * a. The x is obtained
12: from the SNES object (using SNESGetSolution()).
14: Logically Collective on SNES
16: Input Parameters:
17: + snes - the nonlinear solver context
18: . x - the point at which the Jacobian vector products will be performed
19: . jac - the matrix-free Jacobian object
20: . B - either the same as jac or another matrix type (ignored)
21: . flag - not relevent for matrix-free form
22: - dummy - the user context (ignored)
24: Level: developer
26: Warning:
27: If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
28: the x from the SNES object and MatMFFDSetBase() must from that point on be used to
29: change the base vector x.
31: Notes:
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(), MatCreateMFFD(), MATMFFD,
39: MatMFFDSetHHistory(), MatMFFDSetFunctionError(), MatCreateMFFD(), SNESSetJacobian()
41: @*/
42: PetscErrorCode MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
43: {
47: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
48: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
49: return(0);
50: }
52: PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat,MatAssemblyType);
53: PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat,Vec,Vec);
57: /*
58: MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
59: base from the SNES context
61: */
62: PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
63: {
65: MatMFFD j = (MatMFFD)J->data;
66: SNES snes = (SNES)j->ctx;
67: Vec u,f;
70: MatAssemblyEnd_MFFD(J,mt);
72: SNESGetSolution(snes,&u);
73: if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) {
74: SNESGetFunction(snes,&f,NULL,NULL);
75: MatMFFDSetBase_MFFD(J,u,f);
76: } else {
77: /* f value known by SNES is not correct for other differencing function */
78: MatMFFDSetBase_MFFD(J,u,NULL);
79: }
80: return(0);
81: }
83: /*
84: This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
85: uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
86: */
89: PetscErrorCode MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)
90: {
94: MatMFFDSetBase_MFFD(J,U,F);
96: J->ops->assemblyend = MatAssemblyEnd_MFFD;
97: return(0);
98: }
102: /*@
103: MatCreateSNESMF - Creates a matrix-free matrix context for use with
104: a SNES solver. This matrix can be used as the Jacobian argument for
105: the routine SNESSetJacobian(). See MatCreateMFFD() for details on how
106: the finite difference computation is done.
108: Collective on SNES and Vec
110: Input Parameters:
111: . snes - the SNES context
113: Output Parameter:
114: . J - the matrix-free matrix
116: Level: advanced
119: Notes:
120: You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
121: preconditioner matrix
123: If you wish to provide a different function to do differencing on to compute the matrix-free operator than
124: that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call.
126: The difference between this routine and MatCreateMFFD() is that this matrix
127: automatically gets the current base vector from the SNES object and not from an
128: explicit call to MatMFFDSetBase().
130: Warning:
131: If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
132: the x from the SNES object and MatMFFDSetBase() must from that point on be used to
133: change the base vector x.
135: Warning:
136: Using a different function for the differencing will not work if you are using non-linear left preconditioning.
139: .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin()
140: MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(),
141: MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian()
143: @*/
144: PetscErrorCode MatCreateSNESMF(SNES snes,Mat *J)
145: {
147: PetscInt n,N;
148: MatMFFD mf;
151: if (snes->vec_func) {
152: VecGetLocalSize(snes->vec_func,&n);
153: VecGetSize(snes->vec_func,&N);
154: } else if (snes->dm) {
155: Vec tmp;
156: DMGetGlobalVector(snes->dm,&tmp);
157: VecGetLocalSize(tmp,&n);
158: VecGetSize(tmp,&N);
159: DMRestoreGlobalVector(snes->dm,&tmp);
160: } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
161: MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);
162: mf = (MatMFFD)(*J)->data;
163: mf->ctx = snes;
165: if (snes->pc && snes->pcside == PC_LEFT) {
166: MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);
167: } else {
168: MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);
169: }
171: (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;
173: PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);
174: return(0);
175: }