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 */
SNESComputeFunctionCtx(SNES snes,Vec x,Vec f,void * ctx)9 static PetscErrorCode SNESComputeFunctionCtx(SNES snes,Vec x,Vec f,void *ctx)
10 {
11   return SNESComputeFunction(snes,x,f);
12 }
13 
14 /*@C
15     SNESComputeJacobianDefaultColor - Computes the Jacobian using
16     finite differences and coloring to exploit matrix sparsity.
17 
18     Collective on SNES
19 
20     Input Parameters:
21 +   snes - nonlinear solver object
22 .   x1 - location at which to evaluate Jacobian
23 -   ctx - MatFDColoring context or NULL
24 
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)
28 
29     Level: intermediate
30 
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)
37 .  -snes_mf_operator - Use matrix free application of Jacobian
38 -  -snes_mf - Use matrix free Jacobian with not explicit Jacobian represenation
39 
40     Notes:
41         If the coloring is not provided through the context, this will first try to get the
42         coloring from the DM.  If the DM type has no coloring routine, then it will try to
43         get the coloring from the matrix.  This requires that the matrix have nonzero entries
44         precomputed.
45 
46        SNES supports three approaches for computing (approximate) Jacobians: user provided via SNESSetJacobian(), matrix free via SNESSetUseMatrixFree,
47        and computing explictly with finite differences and coloring using MatFDColoring. It is also possible to use automatic differentiation and the MatFDColoring object.
48 
49 
50 .seealso: SNESSetJacobian(), SNESTestJacobian(), SNESComputeJacobianDefault(), SNESSetUseMatrixFree(),
51           MatFDColoringCreate(), MatFDColoringSetFunction()
52 
53 @*/
54 
SNESComputeJacobianDefaultColor(SNES snes,Vec x1,Mat J,Mat B,void * ctx)55 PetscErrorCode  SNESComputeJacobianDefaultColor(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
56 {
57   MatFDColoring  color = (MatFDColoring)ctx;
58   PetscErrorCode ierr;
59   DM             dm;
60   MatColoring    mc;
61   ISColoring     iscoloring;
62   PetscBool      hascolor;
63   PetscBool      solvec,matcolor = PETSC_FALSE;
64 
65   PetscFunctionBegin;
66   if (color) PetscValidHeaderSpecific(color,MAT_FDCOLORING_CLASSID,6);
67   if (!color) {ierr  = PetscObjectQuery((PetscObject)B,"SNESMatFDColoring",(PetscObject*)&color);CHKERRQ(ierr);}
68 
69   if (!color) {
70     ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
71     ierr = DMHasColoring(dm,&hascolor);CHKERRQ(ierr);
72     matcolor = PETSC_FALSE;
73     ierr = PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_fd_color_use_mat",&matcolor,NULL);CHKERRQ(ierr);
74     if (hascolor && !matcolor) {
75       ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&iscoloring);CHKERRQ(ierr);
76       ierr = MatFDColoringCreate(B,iscoloring,&color);CHKERRQ(ierr);
77       ierr = MatFDColoringSetFunction(color,(PetscErrorCode (*)(void))SNESComputeFunctionCtx,NULL);CHKERRQ(ierr);
78       ierr = MatFDColoringSetFromOptions(color);CHKERRQ(ierr);
79       ierr = MatFDColoringSetUp(B,iscoloring,color);CHKERRQ(ierr);
80       ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
81     } else {
82       ierr = MatColoringCreate(B,&mc);CHKERRQ(ierr);
83       ierr = MatColoringSetDistance(mc,2);CHKERRQ(ierr);
84       ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr);
85       ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
86       ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
87       ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
88       ierr = MatFDColoringCreate(B,iscoloring,&color);CHKERRQ(ierr);
89       ierr = MatFDColoringSetFunction(color,(PetscErrorCode (*)(void))SNESComputeFunctionCtx,NULL);CHKERRQ(ierr);
90       ierr = MatFDColoringSetFromOptions(color);CHKERRQ(ierr);
91       ierr = MatFDColoringSetUp(B,iscoloring,color);CHKERRQ(ierr);
92       ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
93     }
94     ierr = PetscObjectCompose((PetscObject)B,"SNESMatFDColoring",(PetscObject)color);CHKERRQ(ierr);
95     ierr = PetscObjectDereference((PetscObject)color);CHKERRQ(ierr);
96   }
97 
98   /* F is only usable if there is no RHS on the SNES and the full solution corresponds to x1 */
99   ierr = VecEqual(x1,snes->vec_sol,&solvec);CHKERRQ(ierr);
100   if (!snes->vec_rhs && solvec) {
101     Vec F;
102     ierr = SNESGetFunction(snes,&F,NULL,NULL);CHKERRQ(ierr);
103     ierr = MatFDColoringSetF(color,F);CHKERRQ(ierr);
104   }
105   ierr = MatFDColoringApply(B,color,x1,snes);CHKERRQ(ierr);
106   if (J != B) {
107     ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
108     ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
109   }
110   PetscFunctionReturn(0);
111 }
112