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: }

 14: /*@C
 15:     SNESComputeJacobianDefaultColor - Computes the Jacobian using
 16:     finite differences and coloring to exploit matrix sparsity.

 18:     Collective on SNES

 20:     Input Parameters:
 21: +   snes - nonlinear solver object
 22: .   x1 - location at which to evaluate Jacobian
 23: -   ctx - MatFDColoring context or NULL

 25:     Output Parameters:
 26: +   J - Jacobian matrix (not altered in this routine)
 27: -   B - newly computed Jacobian matrix to use with preconditioner (generally the same as J)

 29:     Level: intermediate

 31:    Options Database Key:
 32: +  -snes_fd_color_use_mat - use a matrix coloring from the explicit matrix nonzero pattern instead of from the DM providing the matrix
 33: .  -snes_fd_color - Activates SNESComputeJacobianDefaultColor() in SNESSetFromOptions()
 34: .  -mat_fd_coloring_err <err> - Sets <err> (square root of relative error in the function)
 35: .  -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
 36: -  -mat_fd_type - Either wp or ds (see MATMFFD_WP or MATMFFD_DS)

 38:     Notes:
 39:     If the coloring is not provided through the context, this will first try to get the
 40:         coloring from the DM.  If the DM type has no coloring routine, then it will try to
 41:         get the coloring from the matrix.  This requires that the matrix have nonzero entries
 42:         precomputed.  

 44: .seealso: SNESSetJacobian(), SNESTestJacobian(), SNESComputeJacobianDefault()
 45:           MatFDColoringCreate(), MatFDColoringSetFunction()

 47: @*/

 49: PetscErrorCode  SNESComputeJacobianDefaultColor(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
 50: {
 51:   MatFDColoring  color = (MatFDColoring)ctx;
 53:   DM             dm;
 54:   MatColoring    mc;
 55:   ISColoring     iscoloring;
 56:   PetscBool      hascolor;
 57:   PetscBool      solvec,matcolor = PETSC_FALSE;

 61:   if (!color) {PetscObjectQuery((PetscObject)B,"SNESMatFDColoring",(PetscObject*)&color);}

 63:   if (!color) {
 64:     SNESGetDM(snes,&dm);
 65:     DMHasColoring(dm,&hascolor);
 66:     matcolor = PETSC_FALSE;
 67:     PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_fd_color_use_mat",&matcolor,NULL);
 68:     if (hascolor && !matcolor) {
 69:       DMCreateColoring(dm,IS_COLORING_GLOBAL,&iscoloring);
 70:       MatFDColoringCreate(B,iscoloring,&color);
 71:       MatFDColoringSetFunction(color,(PetscErrorCode (*)(void))SNESComputeFunctionCtx,NULL);
 72:       MatFDColoringSetFromOptions(color);
 73:       MatFDColoringSetUp(B,iscoloring,color);
 74:       ISColoringDestroy(&iscoloring);
 75:     } else {
 76:       MatColoringCreate(B,&mc);
 77:       MatColoringSetDistance(mc,2);
 78:       MatColoringSetType(mc,MATCOLORINGSL);
 79:       MatColoringSetFromOptions(mc);
 80:       MatColoringApply(mc,&iscoloring);
 81:       MatColoringDestroy(&mc);
 82:       MatFDColoringCreate(B,iscoloring,&color);
 83:       MatFDColoringSetFunction(color,(PetscErrorCode (*)(void))SNESComputeFunctionCtx,NULL);
 84:       MatFDColoringSetFromOptions(color);
 85:       MatFDColoringSetUp(B,iscoloring,color);
 86:       ISColoringDestroy(&iscoloring);
 87:     }
 88:     PetscObjectCompose((PetscObject)B,"SNESMatFDColoring",(PetscObject)color);
 89:     PetscObjectDereference((PetscObject)color);
 90:   }

 92:   /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
 93:   VecEqual(x1,snes->vec_sol,&solvec);
 94:   if (!snes->vec_rhs && solvec) {
 95:     Vec F;
 96:     SNESGetFunction(snes,&F,NULL,NULL);
 97:     MatFDColoringSetF(color,F);
 98:   }
 99:   MatFDColoringApply(B,color,x1,snes);
100:   if (J != B) {
101:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
102:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
103:   }
104:   return(0);
105: }