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;

 69:   if (!color) PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color);

 71:   if (!color) {
 72:     SNESGetDM(snes, &dm);
 73:     DMHasColoring(dm, &hascolor);
 74:     matcolor = PETSC_FALSE;
 75:     PetscOptionsGetBool(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_fd_color_use_mat", &matcolor, NULL);
 76:     if (hascolor && !matcolor) {
 77:       DMCreateColoring(dm, IS_COLORING_GLOBAL, &iscoloring);
 78:     } else {
 79:       MatColoringCreate(B, &mc);
 80:       MatColoringSetDistance(mc, 2);
 81:       MatColoringSetType(mc, MATCOLORINGSL);
 82:       MatColoringSetFromOptions(mc);
 83:       MatColoringApply(mc, &iscoloring);
 84:       MatColoringDestroy(&mc);
 85:     }
 86:     MatFDColoringCreate(B, iscoloring, &color);
 87:     DMGetDMSNES(dm, &dms);
 88:     if (dms->ops->computemffunction) {
 89:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeMFFunctionCtx, NULL);
 90:     } else {
 91:       MatFDColoringSetFunction(color, (PetscErrorCode(*)(void))SNESComputeFunctionCtx, NULL);
 92:     }
 93:     MatFDColoringSetFromOptions(color);
 94:     MatFDColoringSetUp(B, iscoloring, color);
 95:     ISColoringDestroy(&iscoloring);
 96:     PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)color);
 97:     PetscObjectDereference((PetscObject)color);
 98:   }

100:   /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
101:   VecEqual(x1, snes->vec_sol, &solvec);
102:   if (!snes->vec_rhs && solvec) {
103:     Vec F;
104:     SNESGetFunction(snes, &F, NULL, NULL);
105:     MatFDColoringSetF(color, F);
106:   }
107:   MatFDColoringApply(B, color, x1, snes);
108:   if (J != B) {
109:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
110:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
111:   }
112:   return 0;
113: }