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 { 44 Vec j1a,j2a,x2; 45 PetscInt i,N,start,end,j,value,root,max_funcs = snes->max_funcs; 46 PetscScalar dx,*y,wscale; 47 const PetscScalar *xx; 48 PetscReal amax,epsilon = PETSC_SQRT_MACHINE_EPSILON; 49 PetscReal dx_min = 1.e-16,dx_par = 1.e-1,unorm; 50 MPI_Comm comm; 51 PetscBool assembled,use_wp = PETSC_TRUE,flg; 52 const char *list[2] = {"ds","wp"}; 53 PetscMPIInt size; 54 const PetscInt *ranges; 55 DM dm; 56 DMSNES dms; 57 58 PetscFunctionBegin; 59 snes->max_funcs = PETSC_MAX_INT; 60 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 61 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 62 PetscCall(PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,NULL)); 63 64 PetscCall(PetscObjectGetComm((PetscObject)x1,&comm)); 65 PetscCallMPI(MPI_Comm_size(comm,&size)); 66 PetscCall(MatAssembled(B,&assembled)); 67 if (assembled) PetscCall(MatZeroEntries(B)); 68 if (!snes->nvwork) { 69 if (snes->dm) { 70 PetscCall(DMGetGlobalVector(snes->dm,&j1a)); 71 PetscCall(DMGetGlobalVector(snes->dm,&j2a)); 72 PetscCall(DMGetGlobalVector(snes->dm,&x2)); 73 } else { 74 snes->nvwork = 3; 75 PetscCall(VecDuplicateVecs(x1,snes->nvwork,&snes->vwork)); 76 PetscCall(PetscLogObjectParents(snes,snes->nvwork,snes->vwork)); 77 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 78 } 79 } 80 81 PetscCall(VecGetSize(x1,&N)); 82 PetscCall(VecGetOwnershipRange(x1,&start,&end)); 83 PetscCall(SNESGetDM(snes,&dm)); 84 PetscCall(DMGetDMSNES(dm,&dms)); 85 if (dms->ops->computemffunction) { 86 PetscCall(SNESComputeMFFunction(snes,x1,j1a)); 87 } else { 88 PetscCall(SNESComputeFunction(snes,x1,j1a)); 89 } 90 91 PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES"); 92 PetscCall(PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg)); 93 PetscOptionsEnd(); 94 if (flg && !value) use_wp = PETSC_FALSE; 95 96 if (use_wp) { 97 PetscCall(VecNorm(x1,NORM_2,&unorm)); 98 } 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 PetscCall(VecCopy(x1,x2)); 105 if (i>= start && i<end) { 106 PetscCall(VecGetArrayRead(x1,&xx)); 107 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 108 else dx = xx[i-start]; 109 PetscCall(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 PetscCall(VecSetValues(x2,1,&i,&dx,ADD_VALUES)); 114 } else { 115 wscale = 0.0; 116 } 117 PetscCall(VecAssemblyBegin(x2)); 118 PetscCall(VecAssemblyEnd(x2)); 119 if (dms->ops->computemffunction) { 120 PetscCall(SNESComputeMFFunction(snes,x2,j2a)); 121 } else { 122 PetscCall(SNESComputeFunction(snes,x2,j2a)); 123 } 124 PetscCall(VecAXPY(j2a,-1.0,j1a)); 125 /* Communicate scale=1/dx_i to all processors */ 126 PetscCall(VecGetOwnershipRanges(x1,&ranges)); 127 root = size; 128 for (j=size-1; j>-1; j--) { 129 root--; 130 if (i>=ranges[j]) break; 131 } 132 PetscCallMPI(MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm)); 133 PetscCall(VecScale(j2a,wscale)); 134 PetscCall(VecNorm(j2a,NORM_INFINITY,&amax)); 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) { 138 PetscCall(MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES)); 139 } 140 } 141 PetscCall(VecRestoreArray(j2a,&y)); 142 } 143 if (snes->dm) { 144 PetscCall(DMRestoreGlobalVector(snes->dm,&j1a)); 145 PetscCall(DMRestoreGlobalVector(snes->dm,&j2a)); 146 PetscCall(DMRestoreGlobalVector(snes->dm,&x2)); 147 } 148 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 149 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 150 if (B != J) { 151 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 152 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 153 } 154 snes->max_funcs = max_funcs; 155 snes->nfuncs -= N; 156 PetscFunctionReturn(0); 157 } 158