xref: /petsc/src/snes/mf/snesmfj.c (revision 487a658c8b32ba712a1dc8280daad2fd70c1dcd9)
1 
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>
6 
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()).
11 
12    Logically Collective on SNES
13 
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
18 .   B - either the same as jac or another matrix type (ignored)
19 .   flag - not relevent for matrix-free form
20 -   dummy - the user context (ignored)
21 
22    Level: developer
23 
24    Warning:
25       If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
26     the x from the SNES object and MatMFFDSetBase() must from that point on be used to
27     change the base vector x.
28 
29    Notes:
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.)
35 
36 .seealso: MatMFFDGetH(), MatCreateSNESMF(), MatCreateMFFD(), MATMFFD,
37           MatMFFDSetHHistory(), MatMFFDSetFunctionError(), MatCreateMFFD(), SNESSetJacobian()
38 
39 @*/
40 PetscErrorCode  MatMFFDComputeJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
41 {
42   PetscErrorCode ierr;
43 
44   PetscFunctionBegin;
45   ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
46   ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
47   PetscFunctionReturn(0);
48 }
49 
50 PETSC_EXTERN PetscErrorCode MatAssemblyEnd_MFFD(Mat,MatAssemblyType);
51 PETSC_EXTERN PetscErrorCode MatMFFDSetBase_MFFD(Mat,Vec,Vec);
52 
53 /*@
54     MatSNESMFGetSNES - returns the SNES associated with a matrix created with MatCreateSNESMF()
55 
56     Not collective
57 
58     Input Parameter:
59 .   J - the matrix
60 
61     Output Parameter:
62 .   snes - the SNES object
63 
64 .seealso: MatCreateSNESMF()
65 @*/
66 PetscErrorCode MatSNESMFGetSNES(Mat J,SNES *snes)
67 {
68   MatMFFD        j    = (MatMFFD)J->data;
69 
70   PetscFunctionBegin;
71   *snes = (SNES)j->ctx;
72   PetscFunctionReturn(0);
73 }
74 
75 /*
76    MatAssemblyEnd_SNESMF - Calls MatAssemblyEnd_MFFD() and then sets the
77     base from the SNES context
78 
79 */
80 static PetscErrorCode MatAssemblyEnd_SNESMF(Mat J,MatAssemblyType mt)
81 {
82   PetscErrorCode ierr;
83   MatMFFD        j    = (MatMFFD)J->data;
84   SNES           snes = (SNES)j->ctx;
85   Vec            u,f;
86 
87   PetscFunctionBegin;
88   ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr);
89 
90   ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr);
91   if (j->func == (PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction) {
92     ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);
93     ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr);
94   } else {
95     /* f value known by SNES is not correct for other differencing function */
96     ierr = MatMFFDSetBase_MFFD(J,u,NULL);CHKERRQ(ierr);
97   }
98   PetscFunctionReturn(0);
99 }
100 
101 /*
102    MatAssemblyEnd_SNESMF_UseBase - Calls MatAssemblyEnd_MFFD() and then sets the
103     base from the SNES context. This version will cause the base to be used for differencing
104     even if the func is not SNESComputeFunction. See: MatSNESMFUseBase()
105 
106 
107 */
108 static PetscErrorCode MatAssemblyEnd_SNESMF_UseBase(Mat J,MatAssemblyType mt)
109 {
110   PetscErrorCode ierr;
111   MatMFFD        j    = (MatMFFD)J->data;
112   SNES           snes = (SNES)j->ctx;
113   Vec            u,f;
114 
115   PetscFunctionBegin;
116   ierr = MatAssemblyEnd_MFFD(J,mt);CHKERRQ(ierr);
117 
118   ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr);
119   ierr = SNESGetFunction(snes,&f,NULL,NULL);CHKERRQ(ierr);
120   ierr = MatMFFDSetBase_MFFD(J,u,f);CHKERRQ(ierr);
121   PetscFunctionReturn(0);
122 }
123 
124 /*
125     This routine resets the MatAssemblyEnd() for the MatMFFD created from MatCreateSNESMF() so that it NO longer
126   uses the solution in the SNES object to update the base. See the warning in MatCreateSNESMF().
127 */
128 static PetscErrorCode  MatMFFDSetBase_SNESMF(Mat J,Vec U,Vec F)
129 {
130   PetscErrorCode ierr;
131 
132   PetscFunctionBegin;
133   ierr = MatMFFDSetBase_MFFD(J,U,F);CHKERRQ(ierr);
134   J->ops->assemblyend = MatAssemblyEnd_MFFD;
135   PetscFunctionReturn(0);
136 }
137 
138 static PetscErrorCode  MatSNESMFSetReuseBase_SNESMF(Mat J,PetscBool use)
139 {
140   PetscFunctionBegin;
141   if (use) {
142     J->ops->assemblyend = MatAssemblyEnd_SNESMF_UseBase;
143   } else {
144     J->ops->assemblyend = MatAssemblyEnd_SNESMF;
145   }
146   PetscFunctionReturn(0);
147 }
148 
149 /*@
150     MatSNESMFSetReuseBase - Causes the base vector to be used for differencing even if the function provided to SNESSetFunction() is not the
151                        same as that provided to MatMFFDSetFunction().
152 
153     Logically Collective on Mat
154 
155     Input Parameters:
156 +   J - the MatMFFD matrix
157 -   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is
158           not SNESComputeFunction()
159 
160     Notes:
161     Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with
162            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
163 
164     Developer Notes:
165     This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due
166                      to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration
167                      both functions compute the same mathematical function so the differencing makes sense.
168 
169     Level: advanced
170 
171 .seealso: MatCreateSNESMF(), MatSNESMFGetReuseBase()
172 @*/
173 PetscErrorCode  MatSNESMFSetReuseBase(Mat J,PetscBool use)
174 {
175   PetscErrorCode ierr;
176 
177   PetscFunctionBegin;
178   PetscValidHeaderSpecific(J,MAT_CLASSID,1);
179   ierr = PetscTryMethod(J,"MatSNESMFSetReuseBase_C",(Mat,PetscBool),(J,use));CHKERRQ(ierr);
180   PetscFunctionReturn(0);
181 }
182 
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(0);
189 }
190 
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().
194 
195     Logically Collective on Mat
196 
197     Input Parameter:
198 .   J - the MatMFFD matrix
199 
200     Output Parameter:
201 .   use - if true always reuse the base vector instead of recomputing f(u) even if the function in the MatSNESMF is
202           not SNESComputeFunction()
203 
204     Notes:
205     Care must be taken when using this routine to insure that the function provided to MatMFFDSetFunction(), call it F_MF() is compatible with
206            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
207 
208     Developer Notes:
209     This was provided for the MOOSE team who desired to have a SNESSetFunction() function that could change configurations due
210                      to contacts while the function provided to MatMFFDSetFunction() cannot. Except for the possibility of changing the configuration
211                      both functions compute the same mathematical function so the differencing makes sense.
212 
213     Level: advanced
214 
215 .seealso: MatCreateSNESMF(), MatSNESMFSetReuseBase()
216 @*/
217 PetscErrorCode  MatSNESMFGetReuseBase(Mat J,PetscBool *use)
218 {
219   PetscErrorCode ierr;
220 
221   PetscFunctionBegin;
222   PetscValidHeaderSpecific(J,MAT_CLASSID,1);
223   ierr = PetscUseMethod(J,"MatSNESMFGetReuseBase_C",(Mat,PetscBool*),(J,use));CHKERRQ(ierr);
224   PetscFunctionReturn(0);
225 }
226 
227 /*@
228    MatCreateSNESMF - Creates a matrix-free matrix context for use with
229    a SNES solver.  This matrix can be used as the Jacobian argument for
230    the routine SNESSetJacobian(). See MatCreateMFFD() for details on how
231    the finite difference computation is done.
232 
233    Collective on SNES and Vec
234 
235    Input Parameters:
236 .  snes - the SNES context
237 
238    Output Parameter:
239 .  J - the matrix-free matrix
240 
241    Level: advanced
242 
243 
244    Notes:
245      You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
246      preconditioner matrix
247 
248      If you wish to provide a different function to do differencing on to compute the matrix-free operator than
249      that provided to SNESSetFunction() then call MatMFFDSetFunction() with your function after this call.
250 
251      The difference between this routine and MatCreateMFFD() is that this matrix
252      automatically gets the current base vector from the SNES object and not from an
253      explicit call to MatMFFDSetBase().
254 
255    Warning:
256      If MatMFFDSetBase() is ever called on jac then this routine will NO longer get
257      the x from the SNES object and MatMFFDSetBase() must from that point on be used to
258      change the base vector x.
259 
260    Warning:
261      Using a different function for the differencing will not work if you are using non-linear left preconditioning.
262 
263 
264 .seealso: MatDestroy(), MatMFFDSetFunction(), MatMFFDSetFunctionError(), MatMFFDDSSetUmin()
265           MatMFFDSetHHistory(), MatMFFDResetHHistory(), MatCreateMFFD(),
266           MatMFFDGetH(), MatMFFDRegister(), MatMFFDComputeJacobian(), MatSNESMFSetReuseBase(), MatSNESMFGetReuseBase()
267 
268 @*/
269 PetscErrorCode  MatCreateSNESMF(SNES snes,Mat *J)
270 {
271   PetscErrorCode ierr;
272   PetscInt       n,N;
273   MatMFFD        mf;
274 
275   PetscFunctionBegin;
276   if (snes->vec_func) {
277     ierr = VecGetLocalSize(snes->vec_func,&n);CHKERRQ(ierr);
278     ierr = VecGetSize(snes->vec_func,&N);CHKERRQ(ierr);
279   } else if (snes->dm) {
280     Vec tmp;
281     ierr = DMGetGlobalVector(snes->dm,&tmp);CHKERRQ(ierr);
282     ierr = VecGetLocalSize(tmp,&n);CHKERRQ(ierr);
283     ierr = VecGetSize(tmp,&N);CHKERRQ(ierr);
284     ierr = DMRestoreGlobalVector(snes->dm,&tmp);CHKERRQ(ierr);
285   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() or SNESSetDM() first");
286   ierr = MatCreateMFFD(PetscObjectComm((PetscObject)snes),n,n,N,N,J);CHKERRQ(ierr);
287 
288   mf      = (MatMFFD)(*J)->data;
289   mf->ctx = snes;
290 
291   if (snes->npc && snes->npcside== PC_LEFT) {
292     ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunctionDefaultNPC,snes);CHKERRQ(ierr);
293   } else {
294     ierr = MatMFFDSetFunction(*J,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr);
295   }
296 
297   (*J)->ops->assemblyend = MatAssemblyEnd_SNESMF;
298 
299   ierr = PetscObjectComposeFunction((PetscObject)*J,"MatMFFDSetBase_C",MatMFFDSetBase_SNESMF);CHKERRQ(ierr);
300   ierr = PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFSetReuseBase_C",MatSNESMFSetReuseBase_SNESMF);CHKERRQ(ierr);
301   ierr = PetscObjectComposeFunction((PetscObject)*J,"MatSNESMFGetReuseBase_C",MatSNESMFGetReuseBase_SNESMF);CHKERRQ(ierr);
302   PetscFunctionReturn(0);
303 }
304