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 + x1 - compute Jacobian at this point 12 - ctx - application's function context, as set with SNESSetFunction() 13 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) 17 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) 23 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. 30 31 An alternative routine that uses coloring to exploit matrix sparsity is 32 SNESComputeJacobianDefaultColor(). 33 34 Level: intermediate 35 36 .keywords: SNES, finite differences, Jacobian 37 38 .seealso: SNESSetJacobian(), SNESComputeJacobianDefaultColor(), MatCreateSNESMF() 39 @*/ 40 PetscErrorCode SNESComputeJacobianDefault(SNES snes,Vec x1,Mat J,Mat B,void *ctx) 41 { 42 Vec j1a,j2a,x2; 43 PetscErrorCode ierr; 44 PetscInt i,N,start,end,j,value,root; 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 55 PetscFunctionBegin; 56 /* Since this Jacobian will possibly have "extra" nonzero locations just turn off errors for these locations */ 57 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRQ(ierr); 58 ierr = PetscOptionsGetReal(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_test_err",&epsilon,0);CHKERRQ(ierr); 59 60 ierr = PetscObjectGetComm((PetscObject)x1,&comm);CHKERRQ(ierr); 61 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 62 ierr = MatAssembled(B,&assembled);CHKERRQ(ierr); 63 if (assembled) { 64 ierr = MatZeroEntries(B);CHKERRQ(ierr); 65 } 66 if (!snes->nvwork) { 67 if (snes->dm) { 68 ierr = DMGetGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 69 ierr = DMGetGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 70 ierr = DMGetGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 71 } else { 72 snes->nvwork = 3; 73 ierr = VecDuplicateVecs(x1,snes->nvwork,&snes->vwork);CHKERRQ(ierr); 74 ierr = PetscLogObjectParents(snes,snes->nvwork,snes->vwork);CHKERRQ(ierr); 75 j1a = snes->vwork[0]; j2a = snes->vwork[1]; x2 = snes->vwork[2]; 76 } 77 } 78 79 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 80 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 81 ierr = SNESComputeFunction(snes,x1,j1a);CHKERRQ(ierr); 82 83 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)snes),((PetscObject)snes)->prefix,"Differencing options","SNES");CHKERRQ(ierr); 84 ierr = PetscOptionsEList("-mat_fd_type","Algorithm to compute difference parameter","SNESComputeJacobianDefault",list,2,"wp",&value,&flg);CHKERRQ(ierr); 85 ierr = PetscOptionsEnd();CHKERRQ(ierr); 86 if (flg && !value) use_wp = PETSC_FALSE; 87 88 if (use_wp) { 89 ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr); 90 } 91 /* Compute Jacobian approximation, 1 column at a time. 92 x1 = current iterate, j1a = F(x1) 93 x2 = perturbed iterate, j2a = F(x2) 94 */ 95 for (i=0; i<N; i++) { 96 ierr = VecCopy(x1,x2);CHKERRQ(ierr); 97 if (i>= start && i<end) { 98 ierr = VecGetArrayRead(x1,&xx);CHKERRQ(ierr); 99 if (use_wp) dx = PetscSqrtReal(1.0 + unorm); 100 else dx = xx[i-start]; 101 ierr = VecRestoreArrayRead(x1,&xx);CHKERRQ(ierr); 102 if (PetscAbsScalar(dx) < dx_min) dx = (PetscRealPart(dx) < 0. ? -1. : 1.) * dx_par; 103 dx *= epsilon; 104 wscale = 1.0/dx; 105 ierr = VecSetValues(x2,1,&i,&dx,ADD_VALUES);CHKERRQ(ierr); 106 } else { 107 wscale = 0.0; 108 } 109 ierr = VecAssemblyBegin(x2);CHKERRQ(ierr); 110 ierr = VecAssemblyEnd(x2);CHKERRQ(ierr); 111 ierr = SNESComputeFunction(snes,x2,j2a);CHKERRQ(ierr); 112 ierr = VecAXPY(j2a,-1.0,j1a);CHKERRQ(ierr); 113 /* Communicate scale=1/dx_i to all processors */ 114 ierr = VecGetOwnershipRanges(x1,&ranges);CHKERRQ(ierr); 115 root = size; 116 for (j=size-1; j>-1; j--) { 117 root--; 118 if (i>=ranges[j]) break; 119 } 120 ierr = MPI_Bcast(&wscale,1,MPIU_SCALAR,root,comm);CHKERRQ(ierr); 121 122 ierr = VecScale(j2a,wscale);CHKERRQ(ierr); 123 ierr = VecNorm(j2a,NORM_INFINITY,&amax);CHKERRQ(ierr); amax *= 1.e-14; 124 ierr = VecGetArray(j2a,&y);CHKERRQ(ierr); 125 for (j=start; j<end; j++) { 126 if (PetscAbsScalar(y[j-start]) > amax || j == i) { 127 ierr = MatSetValues(B,1,&j,1,&i,y+j-start,INSERT_VALUES);CHKERRQ(ierr); 128 } 129 } 130 ierr = VecRestoreArray(j2a,&y);CHKERRQ(ierr); 131 } 132 if (snes->dm) { 133 ierr = DMRestoreGlobalVector(snes->dm,&j1a);CHKERRQ(ierr); 134 ierr = DMRestoreGlobalVector(snes->dm,&j2a);CHKERRQ(ierr); 135 ierr = DMRestoreGlobalVector(snes->dm,&x2);CHKERRQ(ierr); 136 } 137 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 138 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 139 if (B != J) { 140 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142 } 143 PetscFunctionReturn(0); 144 } 145 146 147