Actual source code: snesj.c


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

  5: /*@C
  6:    SNESComputeJacobianDefault - Computes the Jacobian using finite differences.

  8:    Collective

 10:    Input Parameters:
 11: +  snes - the `SNES` context
 12: .  x1 - compute Jacobian at this point
 13: -  ctx - application's function context, as set with `SNESSetFunction()`

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

 19:    Options Database Keys:
 20: +  -snes_fd - Activates `SNESComputeJacobianDefault()`
 21: .  -snes_test_err - Square root of function error tolerance, default square root of machine
 22:                     epsilon (1.e-8 in double, 3.e-4 in single)
 23: -  -mat_fd_type - Either wp or ds (see `MATMFFD_WP` or `MATMFFD_DS`)

 25:    Notes:
 26:    This routine is slow and expensive, and is not currently optimized
 27:    to take advantage of sparsity in the problem.  Although
 28:    `SNESComputeJacobianDefault()` is not recommended for general use
 29:    in large-scale applications, It can be useful in checking the
 30:    correctness of a user-provided Jacobian.

 32:    An alternative routine that uses coloring to exploit matrix sparsity is
 33:    `SNESComputeJacobianDefaultColor()`.

 35:    This routine ignores the maximum number of function evaluations set with `SNESSetTolerances()` and the function
 36:    evaluations it performs are not counted in what is returned by of `SNESGetNumberFunctionEvals()`.

 38:    This function can be provided to `SNESSetJacobian()` along with a dense matrix to hold the Jacobian

 40:    Level: intermediate

 42: .seealso: `SNES`, `SNESSetJacobian()`, `SNESSetJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatCreateSNESMF()`
 43: @*/
 44: PetscErrorCode SNESComputeJacobianDefault(SNES snes, Vec x1, Mat J, Mat B, void *ctx)
 45: {
 46:   Vec                j1a, j2a, x2;
 47:   PetscInt           i, N, start, end, j, value, root, max_funcs = snes->max_funcs;
 48:   PetscScalar        dx, *y, wscale;
 49:   const PetscScalar *xx;
 50:   PetscReal          amax, epsilon = PETSC_SQRT_MACHINE_EPSILON;
 51:   PetscReal          dx_min = 1.e-16, dx_par = 1.e-1, unorm;
 52:   MPI_Comm           comm;
 53:   PetscBool          assembled, use_wp = PETSC_TRUE, flg;
 54:   const char        *list[2] = {"ds", "wp"};
 55:   PetscMPIInt        size;
 56:   const PetscInt    *ranges;
 57:   DM                 dm;
 58:   DMSNES             dms;

 60:   snes->max_funcs = PETSC_MAX_INT;
 61:   /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */
 62:   MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 63:   PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_err", &epsilon, NULL);

 65:   PetscObjectGetComm((PetscObject)x1, &comm);
 66:   MPI_Comm_size(comm, &size);
 67:   MatAssembled(B, &assembled);
 68:   if (assembled) MatZeroEntries(B);
 69:   if (!snes->nvwork) {
 70:     if (snes->dm) {
 71:       DMGetGlobalVector(snes->dm, &j1a);
 72:       DMGetGlobalVector(snes->dm, &j2a);
 73:       DMGetGlobalVector(snes->dm, &x2);
 74:     } else {
 75:       snes->nvwork = 3;
 76:       VecDuplicateVecs(x1, snes->nvwork, &snes->vwork);
 77:       j1a = snes->vwork[0];
 78:       j2a = snes->vwork[1];
 79:       x2  = snes->vwork[2];
 80:     }
 81:   }

 83:   VecGetSize(x1, &N);
 84:   VecGetOwnershipRange(x1, &start, &end);
 85:   SNESGetDM(snes, &dm);
 86:   DMGetDMSNES(dm, &dms);
 87:   if (dms->ops->computemffunction) {
 88:     SNESComputeMFFunction(snes, x1, j1a);
 89:   } else {
 90:     SNESComputeFunction(snes, x1, j1a);
 91:   }

 93:   PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing options", "SNES");
 94:   PetscOptionsEList("-mat_fd_type", "Algorithm to compute difference parameter", "SNESComputeJacobianDefault", list, 2, "wp", &value, &flg);
 95:   PetscOptionsEnd();
 96:   if (flg && !value) use_wp = PETSC_FALSE;

 98:   if (use_wp) VecNorm(x1, NORM_2, &unorm);
 99:   /* Compute Jacobian approximation, 1 column at a time.
100:       x1 = current iterate, j1a = F(x1)
101:       x2 = perturbed iterate, j2a = F(x2)
102:    */
103:   for (i = 0; i < N; i++) {
104:     VecCopy(x1, x2);
105:     if (i >= start && i < end) {
106:       VecGetArrayRead(x1, &xx);
107:       if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
108:       else dx = xx[i - start];
109:       VecRestoreArrayRead(x1, &xx);
110:       if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
111:       dx *= epsilon;
112:       wscale = 1.0 / dx;
113:       VecSetValues(x2, 1, &i, &dx, ADD_VALUES);
114:     } else {
115:       wscale = 0.0;
116:     }
117:     VecAssemblyBegin(x2);
118:     VecAssemblyEnd(x2);
119:     if (dms->ops->computemffunction) {
120:       SNESComputeMFFunction(snes, x2, j2a);
121:     } else {
122:       SNESComputeFunction(snes, x2, j2a);
123:     }
124:     VecAXPY(j2a, -1.0, j1a);
125:     /* Communicate scale=1/dx_i to all processors */
126:     VecGetOwnershipRanges(x1, &ranges);
127:     root = size;
128:     for (j = size - 1; j > -1; j--) {
129:       root--;
130:       if (i >= ranges[j]) break;
131:     }
132:     MPI_Bcast(&wscale, 1, MPIU_SCALAR, root, comm);
133:     VecScale(j2a, wscale);
134:     VecNorm(j2a, NORM_INFINITY, &amax);
135:     amax *= 1.e-14;
136:     VecGetArray(j2a, &y);
137:     for (j = start; j < end; j++) {
138:       if (PetscAbsScalar(y[j - start]) > amax || j == i) MatSetValues(B, 1, &j, 1, &i, y + j - start, INSERT_VALUES);
139:     }
140:     VecRestoreArray(j2a, &y);
141:   }
142:   if (snes->dm) {
143:     DMRestoreGlobalVector(snes->dm, &j1a);
144:     DMRestoreGlobalVector(snes->dm, &j2a);
145:     DMRestoreGlobalVector(snes->dm, &x2);
146:   }
147:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
148:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
149:   if (B != J) {
150:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
151:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
152:   }
153:   snes->max_funcs = max_funcs;
154:   snes->nfuncs -= N;
155:   return 0;
156: }