Actual source code: snesj2.c

      1: #include <petsc/private/snesimpl.h>
      2: #include <petscdm.h>

      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(void *snes, Vec x, Vec f, void *ctx)
      9: {
     10:   return SNESComputeFunction((SNES)snes, x, f);
     11: }
     12: static PetscErrorCode SNESComputeMFFunctionCtx(void *snes, Vec x, Vec f, void *ctx)
     13: {
     14:   return SNESComputeMFFunction((SNES)snes, x, f);
     15: }

     17: /*@
     18:   SNESComputeJacobianDefaultColor - Computes the Jacobian using
     19:   finite differences and coloring to exploit matrix sparsity.

     21:   Collective

     23:   Input Parameters:
     24: + snes - nonlinear solver object
     25: . x1   - location at which to evaluate Jacobian
     26: - ctx  - `MatFDColoring` context or `NULL`

     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`)

     32:   Options Database Keys:
     33: + -snes_fd_color_use_mat       - use a matrix coloring from the explicit matrix nonzero pattern instead of from the `DM` providing the matrix
     34: . -snes_fd_color               - Activates `SNESComputeJacobianDefaultColor()` in `SNESSetFromOptions()`
     35: . -mat_fd_coloring_err <err>   - Sets <err> (square root of relative error in the function)
     36: . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude
     37: . -mat_fd_type                 - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)
     38: . -snes_mf_operator            - Use matrix-free application of Jacobian
     39: - -snes_mf                     - Use matrix-free Jacobian with no explicit Jacobian representation

     41:   Notes:

     43:   Level: intermediate

     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:   Developer Note:
     56:   The function has a poorly chosen name since it does not mention the use of finite differences

     58: .seealso: [](ch_snes), `SNES`, `SNESSetJacobian()`, `SNESTestJacobian()`, `SNESComputeJacobianDefault()`, `SNESSetUseMatrixFree()`,
     59:           `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
     60: @*/
     61: PetscErrorCode SNESComputeJacobianDefaultColor(SNES snes, Vec x1, Mat J, Mat B, void *ctx)
     62: {
     63:   MatFDColoring color = (MatFDColoring)ctx;
     64:   DM            dm;
     65:   MatColoring   mc;
     66:   ISColoring    iscoloring;
     67:   PetscBool     hascolor;
     68:   PetscBool     solvec, matcolor = PETSC_FALSE;
     69:   DMSNES        dms;

     71:   PetscFunctionBegin;
     73:   if (!color) PetscCall(PetscObjectQuery((PetscObject)B, "SNESMatFDColoring", (PetscObject *)&color));

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

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

    119: /*@
    120:   SNESPruneJacobianColor - Remove nondiagonal zeros in the Jacobian matrix and update the `MatMFFD` coloring information based on the new nonzero structure

    122:   Collective

    124:   Input Parameters:
    125: + snes - the `SNES` context
    126: . J    - Jacobian matrix (not altered in this routine)
    127: - B    - newly computed Jacobian matrix to use with preconditioner (generally the same as `J`)

    129:   Level: intermediate

    131:   Notes:
    132:   This function improves the `MatMFFD` coloring performance when the Jacobian matrix is overallocated or contains
    133:   many constant zeros entries, which is typically the case when the matrix is generated by a `DM`
    134:   and multiple fields are involved.

    136:   Users need to make sure that the Jacobian matrix is properly filled to reflect the sparsity
    137:   structure. For `MatMFFD` coloring, the values of nonzero entries are not important. So one can
    138:   usually call `SNESComputeJacobian()` with randomized input vectors to generate a dummy Jacobian.
    139:   `SNESComputeJacobian()` should be called before `SNESSolve()` but after `SNESSetUp()`.

    141: .seealso: [](ch_snes), `SNESComputeJacobianDefaultColor()`, `MatEliminateZeros()`, `MatFDColoringCreate()`, `MatFDColoringSetFunction()`
    142: @*/
    143: PetscErrorCode SNESPruneJacobianColor(SNES snes, Mat J, Mat B)
    144: {
    145:   DM            dm;
    146:   DMSNES        dms;
    147:   MatColoring   mc;
    148:   ISColoring    iscoloring;
    149:   MatFDColoring matfdcoloring;

    151:   PetscFunctionBegin;
    152:   /* Generate new coloring after eliminating zeros in the matrix */
    153:   PetscCall(MatEliminateZeros(B, PETSC_TRUE));
    154:   PetscCall(MatColoringCreate(B, &mc));
    155:   PetscCall(MatColoringSetDistance(mc, 2));
    156:   PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
    157:   PetscCall(MatColoringSetFromOptions(mc));
    158:   PetscCall(MatColoringApply(mc, &iscoloring));
    159:   PetscCall(MatColoringDestroy(&mc));
    160:   /* Replace the old coloring with the new one */
    161:   PetscCall(MatFDColoringCreate(B, iscoloring, &matfdcoloring));
    162:   PetscCall(SNESGetDM(snes, &dm));
    163:   PetscCall(DMGetDMSNES(dm, &dms));
    164:   // Comment out the following branch to bypass the coverage test. You can uncomment it when needed.
    165:   //if (dms->ops->computemffunction) {
    166:   //  PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESComputeMFFunctionCtx, NULL));
    167:   //} else {
    168:   PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESComputeFunctionCtx, NULL));
    169:   //}
    170:   PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
    171:   PetscCall(MatFDColoringSetUp(B, iscoloring, matfdcoloring));
    172:   PetscCall(PetscObjectCompose((PetscObject)B, "SNESMatFDColoring", (PetscObject)matfdcoloring));
    173:   PetscCall(PetscObjectDereference((PetscObject)matfdcoloring));
    174:   PetscCall(ISColoringDestroy(&iscoloring));
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }