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

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

  8:    Collective on SNES

 10:    Input Parameters:
 11: +  x1 - compute Jacobian at this point
 12: -  ctx - Section 1.5 Writing Application Codes with PETSc's function context, as set with SNESSetFunction()

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

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

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

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

 34:    Level: intermediate

 36: .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF()
 37: @*/
 38: PetscErrorCode  SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx)
 39: {
 40:   Vec               j1a,j2a,x2;
 41:   PetscErrorCode    ierr;
 42:   PetscInt          i,N,start,end,j,value,root;
 43:   PetscScalar       dx,*y,wscale;
 44:   const PetscScalar *xx;
 45:   PetscReal         amax,epsilon = PETSC_SQRT_MACHINE_EPSILON;
 46:   PetscReal         dx_min = 1.e-16,dx_par = 1.e-1,unorm;
 47:   MPI_Comm          comm;
 48:   PetscBool         assembled,use_wp = PETSC_TRUE,flg;
 49:   const char        *list[2] = {"ds","wp"};
 50:   PetscMPIInt       size;
 51:   const PetscInt    *ranges;

 54:   /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */
 55:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
 56:   PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);

 58:   PetscObjectGetComm((PetscObject)x1,&comm);
 59:   MPI_Comm_size(comm,&size);
 60:   MatAssembled(B,&assembled);
 61:   if (assembled) {
 62:     MatZeroEntries(B);
 63:   }
 64:   if (!snes->nvwork) {
 65:     if (snes->dm) {
 66:       DMGetGlobalVector(snes->dm,&j1a);
 67:       DMGetGlobalVector(snes->dm,&j2a);
 68:       DMGetGlobalVector(snes->dm,&x2);
 69:     } else {
 70:       snes->nvwork = 3;
 71:       VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);
 72:       PetscLogObjectParents(snes,snes->nvwork,snes->vwork);
 73:       j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2];
 74:     }
 75:   }

 77:   VecGetSize(x1,&N);
 78:   VecGetOwnershipRange(x1,&start,&end);
 79:   SNESComputeFunction(snes,x1,j1a);

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

 86:   if (use_wp) {
 87:     VecNorm(x1,NORM_2,&unorm);
 88:   }
 89:   /* Compute Jacobian approximation, 1 column at a time.
 90:       x1 = current iterate, j1a = F(x1)
 91:       x2 = perturbed iterate, j2a = F(x2)
 92:    */
 93:   for (i=0; i<N; i++) {
 94:     VecCopy(x1,x2);
 95:     if (i>= start && i<end) {
 96:       VecGetArrayRead(x1,&xx);
 97:       if (use_wp) dx = PetscSqrtReal(1.0 + unorm);
 98:       else        dx = xx[i-start];
 99:       VecRestoreArrayRead(x1,&xx);
100:       if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par;
101:       dx    *= epsilon;
102:       wscale = 1.0/dx;
103:       VecSetValues(x2,1,&i,&dx,ADD_VALUES);
104:     } else {
105:       wscale = 0.0;
106:     }
107:     VecAssemblyBegin(x2);
108:     VecAssemblyEnd(x2);
109:     SNESComputeFunction(snes,x2,j2a);
110:     VecAXPY(j2a,-1.0,j1a);
111:     /* Communicate scale=1/dx_i to all processors */
112:     VecGetOwnershipRanges(x1,&ranges);
113:     root = size;
114:     for (j=size-1; j>-1; j--) {
115:       root--;
116:       if (i>=ranges[j]) break;
117:     }
118:     MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);

120:     VecScale(j2a,wscale);
121:     VecNorm(j2a,NORM_INFINITY,&amax); amax *= 1.e-14;
122:     VecGetArray(j2a,&y);
123:     for (j=start; j<end; j++) {
124:       if (PetscAbsScalar(y[j-start]) > amax || j == i) {
125:         MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);
126:       }
127:     }
128:     VecRestoreArray(j2a,&y);
129:   }
130:   if (snes->dm) {
131:     DMRestoreGlobalVector(snes->dm,&j1a);
132:     DMRestoreGlobalVector(snes->dm,&j2a);
133:     DMRestoreGlobalVector(snes->dm,&x2);
134:   }
135:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
136:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
137:   if (B != J) {
138:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
139:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
140:   }
141:   return(0);
142: }