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 */ 8 static PetscErrorCode SNESComputeFunctionCtx(void *snes, Vec x, Vec f, void *ctx) 9 { 10 return SNESComputeFunction((SNES)snes, x, f); 11 } 12 static PetscErrorCode SNESComputeMFFunctionCtx(void *snes, Vec x, Vec f, void *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 @*/ 61 PetscErrorCode SNESComputeJacobianDefaultColor(SNES snes, Vec x1, Mat J, Mat B, void *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 @*/ 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