1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petscdm.h> /*I "petscdm.h" I*/
3
4 /*
5 MatFDColoringSetFunction() takes a function with four arguments, we want to use SNESComputeFunction()
6 since it logs function computation information.
7 */
SNESComputeFunctionCtx(void * snes,Vec x,Vec f,PetscCtx ctx)8 static PetscErrorCode SNESComputeFunctionCtx(void *snes, Vec x, Vec f, PetscCtx ctx)
9 {
10 return SNESComputeFunction((SNES)snes, x, f);
11 }
SNESComputeMFFunctionCtx(void * snes,Vec x,Vec f,PetscCtx ctx)12 static PetscErrorCode SNESComputeMFFunctionCtx(void *snes, Vec x, Vec f, PetscCtx ctx)
13 {
14 return SNESComputeMFFunction((SNES)snes, x, f);
15 }
16
17 /*@
18 SNESComputeJacobianDefaultColor - Computes the Jacobian using
19 finite differences and coloring to exploit matrix sparsity.
20
21 Collective
22
23 Input Parameters:
24 + snes - nonlinear solver object
25 . x1 - location at which to evaluate Jacobian
26 - ctx - `MatFDColoring` context or `NULL`
27
28 Output Parameters:
29 + J - Jacobian matrix (not altered in this routine)
30 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
31
32 Options Database Keys:
33 + -snes_fd_color_use_mat - use a matrix coloring from the explicit matrix nonzero pattern instead of from the `DM` providing the matrix
34 . -snes_fd_color - Activates `SNESComputeJacobianDefaultColor()` in `SNESSetFromOptions()`
35 . -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
36 . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
37 . -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
38 . -snes_mf_operator - Use matrix-free application of Jacobian
39 - -snes_mf - Use matrix-free Jacobian with no explicit Jacobian representation
40
41 Notes:
42
43 Level: intermediate
44
45 If the coloring is not provided through the context, this will first try to get the
46 coloring from the `DM`. If the `DM` has no coloring routine, then it will try to
47 get the coloring from the matrix. This requires that the matrix have its nonzero locations already provided.
48
49 `SNES` supports three approaches for computing (approximate) Jacobians: user provided via `SNESSetJacobian()`, matrix-free via `SNESSetUseMatrixFree()`,
50 and computing explicitly with finite differences and coloring using `MatFDColoring`. It is also possible to use automatic differentiation
51 and the `MatFDColoring` object, see src/ts/tutorials/autodiff/ex16adj_tl.cxx
52
53 This function can be provided to `SNESSetJacobian()` along with an appropriate sparse matrix to hold the Jacobian
54
55 Developer Note:
56 The function has a poorly chosen name since it does not mention the use of finite differences
57
58 .seealso: [](ch_snes), `SNES`, `SNESSetJacobian()`, `SNESTestJacobian()`, `SNESComputeJacobianDefault()`, `SNESSetUseMatrixFree()`,
59 `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
60 @*/
SNESComputeJacobianDefaultColor(SNES snes,Vec x1,Mat J,Mat B,PetscCtx ctx)61 PetscErrorCode SNESComputeJacobianDefaultColor(SNES snes, Vec x1, Mat J, Mat B, PetscCtx ctx)
62 {
63 MatFDColoring color = (MatFDColoring)ctx;
64 DM dm;
65 MatColoring mc;
66 ISColoring iscoloring;
67 PetscBool hascolor;
68 PetscBool solvec, matcolor = PETSC_FALSE;
69 DMSNES dms;
70
71 PetscFunctionBegin;
72 if (color) PetscValidHeaderSpecific(color, MAT_FDCOLORING_CLASSID, 5);
73 if (!color) PetscCall(PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color));
74
75 if (!color) {
76 PetscCall(SNESGetDM(snes, &dm));
77 PetscCall(DMHasColoring(dm, &hascolor));
78 matcolor = PETSC_FALSE;
79 PetscCall(PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fd_color_use_mat", &matcolor, NULL));
80 if (hascolor && !matcolor) {
81 PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
82 } else {
83 PetscCall(MatColoringCreate(B, &mc));
84 PetscCall(MatColoringSetDistance(mc, 2));
85 PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
86 PetscCall(MatColoringSetFromOptions(mc));
87 PetscCall(MatColoringApply(mc, &iscoloring));
88 PetscCall(MatColoringDestroy(&mc));
89 }
90 PetscCall(MatFDColoringCreate(B, iscoloring, &color));
91 PetscCall(DMGetDMSNES(dm, &dms));
92 if (dms->ops->computemffunction) {
93 PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESComputeMFFunctionCtx, NULL));
94 } else {
95 PetscCall(MatFDColoringSetFunction(color, (MatFDColoringFn *)SNESComputeFunctionCtx, NULL));
96 }
97 PetscCall(MatFDColoringSetFromOptions(color));
98 PetscCall(MatFDColoringSetUp(B, iscoloring, color));
99 PetscCall(ISColoringDestroy(&iscoloring));
100 PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)color));
101 PetscCall(PetscObjectDereference((PetscObject)color));
102 }
103
104 /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
105 PetscCall(VecEqual(x1, snes->vec_sol, &solvec));
106 if (!snes->vec_rhs && solvec) {
107 Vec F;
108 PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
109 PetscCall(MatFDColoringSetF(color, F));
110 }
111 PetscCall(MatFDColoringApply(B, color, x1, snes));
112 if (J != B) {
113 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
114 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
115 }
116 PetscFunctionReturn(PETSC_SUCCESS);
117 }
118
119 /*@
120 SNESPruneJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information based on the new nonzero structure
121
122 Collective
123
124 Input Parameters:
125 + snes - the `SNES` context
126 . J - Jacobian matrix (not altered in this routine)
127 - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
128
129 Level: intermediate
130
131 Notes:
132 This function improves the `MatMFFD` coloring performance when the Jacobian matrix is overallocated or contains
133 many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
134 and multiple fields are involved.
135
136 Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
137 structure. For `MatMFFD` coloring, the values of nonzero entries are not important. So one can
138 usually call `SNESComputeJacobian()` with randomized input vectors to generate a dummy Jacobian.
139 `SNESComputeJacobian()` should be called before `SNESSolve()` but after `SNESSetUp()`.
140
141 .seealso: [](ch_snes), `SNESComputeJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
142 @*/
SNESPruneJacobianColor(SNES snes,Mat J,Mat B)143 PetscErrorCode SNESPruneJacobianColor(SNES snes, Mat J, Mat B)
144 {
145 DM dm;
146 DMSNES dms;
147 MatColoring mc;
148 ISColoring iscoloring;
149 MatFDColoring matfdcoloring;
150
151 PetscFunctionBegin;
152 /* Generate new coloring after eliminating zeros in the matrix */
153 PetscCall(MatEliminateZeros(B, PETSC_TRUE));
154 PetscCall(MatColoringCreate(B, &mc));
155 PetscCall(MatColoringSetDistance(mc, 2));
156 PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
157 PetscCall(MatColoringSetFromOptions(mc));
158 PetscCall(MatColoringApply(mc, &iscoloring));
159 PetscCall(MatColoringDestroy(&mc));
160 /* Replace the old coloring with the new one */
161 PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
162 PetscCall(SNESGetDM(snes, &dm));
163 PetscCall(DMGetDMSNES(dm, &dms));
164 // Comment out the following branch to bypass the coverage test. You can uncomment it when needed.
165 //if (dms->ops->computemffunction) {
166 // PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESComputeMFFunctionCtx, NULL));
167 //} else {
168 PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESComputeFunctionCtx, NULL));
169 //}
170 PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
171 PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
172 PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)matfdcoloring));
173 PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
174 PetscCall(ISColoringDestroy(&iscoloring));
175 PetscFunctionReturn(PETSC_SUCCESS);
176 }
177