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