1 2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdm.h> 4 5 /*@C 6 SNESComputeJacobianDefault - Computes the Jacobian using finite differences. 7 8 Collective on SNES 9 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() 14 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) 18 19 Options Database Key: 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) 24 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. 31 32 An alternative routine that uses coloring to exploit matrix sparsity is 33 SNESComputeJacobianDefaultColor(). 34 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(). 37 38 Level: intermediate 39 40 .seealso: `SNESSetJacobian()`, `SNESComputeJacobianDefaultColor()`, `MatCreateSNESMF()` 41 @*/ 42 PetscErrorCode SNESComputeJacobianDefault(SNES snes, Vec x1, Mat J, Mat B, void *ctx) { 43 Vec j1a, j2a, x2; 44 PetscInt i, N, start, end, j, value, root, max_funcs = snes->max_funcs; 45 PetscScalar dx, *y, wscale; 46 const PetscScalar *xx; 47 PetscReal amax, epsilon = PETSC_SQRT_MACHINE_EPSILON; 48 PetscReal dx_min = 1.e-16, dx_par = 1.e-1, unorm; 49 MPI_Comm comm; 50 PetscBool assembled, use_wp = PETSC_TRUE, flg; 51 const char *list[2] = {"ds", "wp"}; 52 PetscMPIInt size; 53 const PetscInt *ranges; 54 DM dm; 55 DMSNES dms; 56 57 PetscFunctionBegin; 58 snes->max_funcs = PETSC_MAX_INT; 59 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 60 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 61 PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_test_err", &epsilon, NULL)); 62 63 PetscCall(PetscObjectGetComm((PetscObject)x1, &comm)); 64 PetscCallMPI(MPI_Comm_size(comm, &size)); 65 PetscCall(MatAssembled(B, &assembled)); 66 if (assembled) PetscCall(MatZeroEntries(B)); 67 if (!snes->nvwork) { 68 if (snes->dm) { 69 PetscCall(DMGetGlobalVector(snes->dm, &j1a)); 70 PetscCall(DMGetGlobalVector(snes->dm, &j2a)); 71 PetscCall(DMGetGlobalVector(snes->dm, &x2)); 72 } else { 73 snes->nvwork = 3; 74 PetscCall(VecDuplicateVecs(x1, snes->nvwork, &snes->vwork)); 75 PetscCall(PetscLogObjectParents(snes, snes->nvwork, snes->vwork)); 76 j1a = snes->vwork[0]; 77 j2a = snes->vwork[1]; 78 x2 = snes->vwork[2]; 79 } 80 } 81 82 PetscCall(VecGetSize(x1, &N)); 83 PetscCall(VecGetOwnershipRange(x1, &start, &end)); 84 PetscCall(SNESGetDM(snes, &dm)); 85 PetscCall(DMGetDMSNES(dm, &dms)); 86 if (dms->ops->computemffunction) { 87 PetscCall(SNESComputeMFFunction(snes, x1, j1a)); 88 } else { 89 PetscCall(SNESComputeFunction(snes, x1, j1a)); 90 } 91 92 PetscOptionsBegin(PetscObjectComm((PetscObject)snes), ((PetscObject)snes)->prefix, "Differencing options", "SNES"); 93 PetscCall(PetscOptionsEList("-mat_fd_type", "Algorithm to compute difference parameter", "SNESComputeJacobianDefault", list, 2, "wp", &value, &flg)); 94 PetscOptionsEnd(); 95 if (flg && !value) use_wp = PETSC_FALSE; 96 97 if (use_wp) PetscCall(VecNorm(x1, NORM_2, &unorm)); 98 /* Compute Jacobian approximation, 1 column at a time. 99 x1 = current iterate, j1a = F(x1) 100 x2 = perturbed iterate, j2a = F(x2) 101 */ 102 for (i = 0; i < N; i++) { 103 PetscCall(VecCopy(x1, x2)); 104 if (i >= start && i < end) { 105 PetscCall(VecGetArrayRead(x1, &xx)); 106 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 107 else dx = xx[i - start]; 108 PetscCall(VecRestoreArrayRead(x1, &xx)); 109 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 110 dx *= epsilon; 111 wscale = 1.0 / dx; 112 PetscCall(VecSetValues(x2, 1, &i, &dx, ADD_VALUES)); 113 } else { 114 wscale = 0.0; 115 } 116 PetscCall(VecAssemblyBegin(x2)); 117 PetscCall(VecAssemblyEnd(x2)); 118 if (dms->ops->computemffunction) { 119 PetscCall(SNESComputeMFFunction(snes, x2, j2a)); 120 } else { 121 PetscCall(SNESComputeFunction(snes, x2, j2a)); 122 } 123 PetscCall(VecAXPY(j2a, -1.0, j1a)); 124 /* Communicate scale=1/dx_i to all processors */ 125 PetscCall(VecGetOwnershipRanges(x1, &ranges)); 126 root = size; 127 for (j = size - 1; j > -1; j--) { 128 root--; 129 if (i >= ranges[j]) break; 130 } 131 PetscCallMPI(MPI_Bcast(&wscale, 1, MPIU_SCALAR, root, comm)); 132 PetscCall(VecScale(j2a, wscale)); 133 PetscCall(VecNorm(j2a, NORM_INFINITY, &amax)); 134 amax *= 1.e-14; 135 PetscCall(VecGetArray(j2a, &y)); 136 for (j = start; j < end; j++) { 137 if (PetscAbsScalar(y[j - start]) > amax || j == i) PetscCall(MatSetValues(B, 1, &j, 1, &i, y + j - start, INSERT_VALUES)); 138 } 139 PetscCall(VecRestoreArray(j2a, &y)); 140 } 141 if (snes->dm) { 142 PetscCall(DMRestoreGlobalVector(snes->dm, &j1a)); 143 PetscCall(DMRestoreGlobalVector(snes->dm, &j2a)); 144 PetscCall(DMRestoreGlobalVector(snes->dm, &x2)); 145 } 146 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 147 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 148 if (B != J) { 149 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 150 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 151 } 152 snes->max_funcs = max_funcs; 153 snes->nfuncs -= N; 154 PetscFunctionReturn(0); 155 } 156