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 { 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 } 17 18 /*@C 19 SNESComputeJacobianDefaultColor - Computes the Jacobian using 20 finite differences and coloring to exploit matrix sparsity. 21 22 Collective 23 24 Input Parameters: 25 + snes - nonlinear solver object 26 . x1 - location at which to evaluate Jacobian 27 - ctx - `MatFDColoring` context or `NULL` 28 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`) 32 33 Level: intermediate 34 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 43 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. 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 .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; 67 68 PetscFunctionBegin; 69 if (color) PetscValidHeaderSpecific(color, MAT_FDCOLORING_CLASSID, 5); 70 if (!color) PetscCall(PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color)); 71 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 } 100 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 } 115 116 /*@C 117 SNESPruneJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information. 118 119 Collective 120 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`) 125 126 Level: intermediate 127 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. 132 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()`. 137 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; 147 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 } 174