xref: /petsc/src/snes/interface/snesj2.c (revision 6a4a1270853b5b1995c94de98f44d28bec57470a)
1 
2 #include <petsc/private/snesimpl.h> /*I  "petscsnes.h"  I*/
3 #include <petscdm.h>                /*I  "petscdm.h"    I*/
4 
5 /*
6    MatFDColoringSetFunction() takes a function with four arguments, we want to use SNESComputeFunction()
7    since it logs function computation information.
8 */
9 static PetscErrorCode SNESComputeFunctionCtx(SNES snes, Vec x, Vec f, void *ctx) {
10   return SNESComputeFunction(snes, x, f);
11 }
12 static PetscErrorCode SNESComputeMFFunctionCtx(SNES snes, Vec x, Vec f, void *ctx) {
13   return SNESComputeMFFunction(snes, x, f);
14 }
15 
16 /*@C
17     SNESComputeJacobianDefaultColor - Computes the Jacobian using
18     finite differences and coloring to exploit matrix sparsity.
19 
20     Collective on snes
21 
22     Input Parameters:
23 +   snes - nonlinear solver object
24 .   x1 - location at which to evaluate Jacobian
25 -   ctx - MatFDColoring context or NULL
26 
27     Output Parameters:
28 +   J - Jacobian matrix (not altered in this routine)
29 -   B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)
30 
31     Level: intermediate
32 
33    Options Database Keys:
34 +  -snes_fd_color_use_mat - use a matrix coloring from the explicit matrix nonzero pattern instead of from the DM providing the matrix
35 .  -snes_fd_color - Activates `SNESComputeJacobianDefaultColor()` in `SNESSetFromOptions()`
36 .  -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
37 .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
38 .  -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
39 .  -snes_mf_operator - Use matrix free application of Jacobian
40 -  -snes_mf - Use matrix free Jacobian with no explicit Jacobian representation
41 
42     Notes:
43     If the coloring is not provided through the context, this will first try to get the
44     coloring from the `DM`.  If the `DM` has no coloring routine, then it will try to
45     get the coloring from the matrix.  This requires that the matrix have its nonzero locations already provided.
46 
47     `SNES` supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix-free via `SNESSetUseMatrixFree()`,
48     and computing explicitly with finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation
49     and the `MatFDColoring` object, see src/ts/tutorials/autodiff/ex16adj_tl.cxx
50 
51    This function can be provided to `SNESSetJacobian()` along with an appropriate sparse matrix to hold the Jacobian
52 
53 .seealso: `SNES`, `SNESSetJacobian()`, `SNESTestJacobian()`, `SNESComputeJacobianDefault()`, `SNESSetUseMatrixFree()`,
54           `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
55 @*/
56 PetscErrorCode SNESComputeJacobianDefaultColor(SNES snes, Vec x1, Mat J, Mat B, void *ctx) {
57   MatFDColoring color = (MatFDColoring)ctx;
58   DM            dm;
59   MatColoring   mc;
60   ISColoring    iscoloring;
61   PetscBool     hascolor;
62   PetscBool     solvec, matcolor = PETSC_FALSE;
63   DMSNES        dms;
64 
65   PetscFunctionBegin;
66   if (color) PetscValidHeaderSpecific(color, MAT_FDCOLORING_CLASSID, 5);
67   if (!color) PetscCall(PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color));
68 
69   if (!color) {
70     PetscCall(SNESGetDM(snes, &dm));
71     PetscCall(DMHasColoring(dm, &hascolor));
72     matcolor = PETSC_FALSE;
73     PetscCall(PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fd_color_use_mat", &matcolor, NULL));
74     if (hascolor && !matcolor) {
75       PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
76     } else {
77       PetscCall(MatColoringCreate(B, &mc));
78       PetscCall(MatColoringSetDistance(mc, 2));
79       PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
80       PetscCall(MatColoringSetFromOptions(mc));
81       PetscCall(MatColoringApply(mc, &iscoloring));
82       PetscCall(MatColoringDestroy(&mc));
83     }
84     PetscCall(MatFDColoringCreate(B, iscoloring, &color));
85     PetscCall(DMGetDMSNES(dm, &dms));
86     if (dms->ops->computemffunction) {
87       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeMFFunctionCtx, NULL));
88     } else {
89       PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeFunctionCtx, NULL));
90     }
91     PetscCall(MatFDColoringSetFromOptions(color));
92     PetscCall(MatFDColoringSetUp(B, iscoloring, color));
93     PetscCall(ISColoringDestroy(&iscoloring));
94     PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)color));
95     PetscCall(PetscObjectDereference((PetscObject)color));
96   }
97 
98   /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
99   PetscCall(VecEqual(x1, snes->vec_sol, &solvec));
100   if (!snes->vec_rhs && solvec) {
101     Vec F;
102     PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
103     PetscCall(MatFDColoringSetF(color, F));
104   }
105   PetscCall(MatFDColoringApply(B, color, x1, snes));
106   if (J != B) {
107     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
108     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
109   }
110   PetscFunctionReturn(0);
111 }
112