Actual source code: snesj2.c
2: #include <petsc/private/snesimpl.h>
3: #include <petscdm.h>
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: {
11: return SNESComputeFunction(snes, x, f);
12: }
13: static PetscErrorCode SNESComputeMFFunctionCtx(SNES snes, Vec x, Vec f, void *ctx)
14: {
15: return SNESComputeMFFunction(snes, x, f);
16: }
18: /*@C
19: SNESComputeJacobianDefaultColor - Computes the Jacobian using
20: finite differences and coloring to exploit matrix sparsity.
22: Collective
24: Input Parameters:
25: + snes - nonlinear solver object
26: . x1 - location at which to evaluate Jacobian
27: - ctx - `MatFDColoring` context or `NULL`
29: Output Parameters:
30: + J - Jacobian matrix (not altered in this routine)
31: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
33: Level: intermediate
35: Options Database Keys:
36: + -snes_fd_color_use_mat - use a matrix coloring from the explicit matrix nonzero pattern instead of from the `DM` providing the matrix
37: . -snes_fd_color - Activates `SNESComputeJacobianDefaultColor()` in `SNESSetFromOptions()`
38: . -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
39: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
40: . -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
41: . -snes_mf_operator - Use matrix free application of Jacobian
42: - -snes_mf - Use matrix free Jacobian with no explicit Jacobian representation
44: Notes:
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.
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
53: This function can be provided to `SNESSetJacobian()` along with an appropriate sparse matrix to hold the Jacobian
55: .seealso: `SNES`, `SNESSetJacobian()`, `SNESTestJacobian()`, `SNESComputeJacobianDefault()`, `SNESSetUseMatrixFree()`,
56: `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
57: @*/
58: PetscErrorCode SNESComputeJacobianDefaultColor(SNES snes, Vec x1, Mat J, Mat B, void *ctx)
59: {
60: MatFDColoring color = (MatFDColoring)ctx;
61: DM dm;
62: MatColoring mc;
63: ISColoring iscoloring;
64: PetscBool hascolor;
65: PetscBool solvec, matcolor = PETSC_FALSE;
66: DMSNES dms;
68: PetscFunctionBegin;
70: if (!color) PetscCall(PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color));
72: if (!color) {
73: PetscCall(SNESGetDM(snes, &dm));
74: PetscCall(DMHasColoring(dm, &hascolor));
75: matcolor = PETSC_FALSE;
76: PetscCall(PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fd_color_use_mat", &matcolor, NULL));
77: if (hascolor && !matcolor) {
78: PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring));
79: } else {
80: PetscCall(MatColoringCreate(B, &mc));
81: PetscCall(MatColoringSetDistance(mc, 2));
82: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
83: PetscCall(MatColoringSetFromOptions(mc));
84: PetscCall(MatColoringApply(mc, &iscoloring));
85: PetscCall(MatColoringDestroy(&mc));
86: }
87: PetscCall(MatFDColoringCreate(B, iscoloring, &color));
88: PetscCall(DMGetDMSNES(dm, &dms));
89: if (dms->ops->computemffunction) {
90: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeMFFunctionCtx, NULL));
91: } else {
92: PetscCall(MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeFunctionCtx, NULL));
93: }
94: PetscCall(MatFDColoringSetFromOptions(color));
95: PetscCall(MatFDColoringSetUp(B, iscoloring, color));
96: PetscCall(ISColoringDestroy(&iscoloring));
97: PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)color));
98: PetscCall(PetscObjectDereference((PetscObject)color));
99: }
101: /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
102: PetscCall(VecEqual(x1, snes->vec_sol, &solvec));
103: if (!snes->vec_rhs && solvec) {
104: Vec F;
105: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
106: PetscCall(MatFDColoringSetF(color, F));
107: }
108: PetscCall(MatFDColoringApply(B, color, x1, snes));
109: if (J != B) {
110: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
111: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
112: }
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@C
117: SNESPruneJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information.
119: Collective
121: Input Parameters:
122: + snes - the `SNES` context
123: . J - Jacobian matrix (not altered in this routine)
124: - B - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)
126: Level: intermediate
128: Notes:
129: This function improves the `MatMFFD` coloring performance when the Jacobian matrix is overallocated or contains
130: many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
131: and multiple fields are involved.
133: Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
134: structure. For `MatMFFD` coloring, the values of nonzero entries are not important. So one can
135: usually call `SNESComputeJacobian()` with randomized input vectors to generate a dummy Jacobian.
136: `SNESComputeJacobian()` should be called before `SNESSolve()` but after `SNESSetUp()`.
138: .seealso: `SNESComputeJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
139: @*/
140: PetscErrorCode SNESPruneJacobianColor(SNES snes, Mat J, Mat B)
141: {
142: DM dm;
143: DMSNES dms;
144: MatColoring mc;
145: ISColoring iscoloring;
146: MatFDColoring matfdcoloring;
148: PetscFunctionBegin;
149: /* Generate new coloring after eliminating zeros in the matrix */
150: PetscCall(MatEliminateZeros(B));
151: PetscCall(MatColoringCreate(B, &mc));
152: PetscCall(MatColoringSetDistance(mc, 2));
153: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
154: PetscCall(MatColoringSetFromOptions(mc));
155: PetscCall(MatColoringApply(mc, &iscoloring));
156: PetscCall(MatColoringDestroy(&mc));
157: /* Replace the old coloring with the new one */
158: PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
159: PetscCall(SNESGetDM(snes, &dm));
160: PetscCall(DMGetDMSNES(dm, &dms));
161: // Comment out the following branch to bypass the coverage test. You can uncomment it when needed.
162: //if (dms->ops->computemffunction) {
163: // PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESComputeMFFunctionCtx, NULL));
164: //} else {
165: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESComputeFunctionCtx, NULL));
166: //}
167: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
168: PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
169: PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)matfdcoloring));
170: PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
171: PetscCall(ISColoringDestroy(&iscoloring));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }