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 on snes 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(0); 114 } 115 116 /*@C 117 SNESPruneJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the MFFD coloring. 118 119 Collective on SNES 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 MFFD coloring 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 MFFD 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(0); 173 } 174