1 2 static const char help[] = "Tries to solve u`` + u^{2} = f for an easy case and an impossible case.\n\n"; 3 4 /* 5 This example was contributed by Peter Graf to show how SNES fails when given a nonlinear problem with no solution. 6 7 Run with -n 14 to see fail to converge and -n 15 to see convergence 8 9 The option -second_order uses a different discretization of the Neumann boundary condition and always converges 10 11 */ 12 13 #include <petscsnes.h> 14 15 PetscBool second_order = PETSC_FALSE; 16 #define X0DOT -2.0 17 #define X1 5.0 18 #define KPOW 2.0 19 const PetscScalar sperturb = 1.1; 20 21 /* 22 User-defined routines 23 */ 24 PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 25 PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 26 27 int main(int argc,char **argv) 28 { 29 SNES snes; /* SNES context */ 30 Vec x,r,F; /* vectors */ 31 Mat J; /* Jacobian */ 32 PetscInt it,n = 11,i; 33 PetscReal h,xp = 0.0; 34 PetscScalar v; 35 const PetscScalar a = X0DOT; 36 const PetscScalar b = X1; 37 const PetscScalar k = KPOW; 38 PetscScalar v2; 39 PetscScalar *xx; 40 41 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 42 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 43 PetscCall(PetscOptionsGetBool(NULL,NULL,"-second_order",&second_order,NULL)); 44 h = 1.0/(n-1); 45 46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 Create nonlinear solver context 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 50 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create vector data structures; set function evaluation routine 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 56 PetscCall(VecCreate(PETSC_COMM_SELF,&x)); 57 PetscCall(VecSetSizes(x,PETSC_DECIDE,n)); 58 PetscCall(VecSetFromOptions(x)); 59 PetscCall(VecDuplicate(x,&r)); 60 PetscCall(VecDuplicate(x,&F)); 61 62 PetscCall(SNESSetFunction(snes,r,FormFunction,(void*)F)); 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Create matrix data structures; set Jacobian evaluation routine 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 67 68 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,&J)); 69 70 /* 71 Note that in this case we create separate matrices for the Jacobian 72 and preconditioner matrix. Both of these are computed in the 73 routine FormJacobian() 74 */ 75 /* PetscCall(SNESSetJacobian(snes,NULL,JPrec,FormJacobian,0)); */ 76 PetscCall(SNESSetJacobian(snes,J,J,FormJacobian,0)); 77 /* PetscCall(SNESSetJacobian(snes,J,JPrec,FormJacobian,0)); */ 78 79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80 Customize nonlinear solver; set runtime options 81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 82 83 PetscCall(SNESSetFromOptions(snes)); 84 85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 Initialize application: 87 Store right-hand-side of PDE and exact solution 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 90 /* set right hand side and initial guess to be exact solution of continuim problem */ 91 #define SQR(x) ((x)*(x)) 92 xp = 0.0; 93 for (i=0; i<n; i++) 94 { 95 v = k*(k-1.)*(b-a)*PetscPowScalar(xp,k-2.) + SQR(a*xp) + SQR(b-a)*PetscPowScalar(xp,2.*k) + 2.*a*(b-a)*PetscPowScalar(xp,k+1.); 96 PetscCall(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 97 v2 = a*xp + (b-a)*PetscPowScalar(xp,k); 98 PetscCall(VecSetValues(x,1,&i,&v2,INSERT_VALUES)); 99 xp += h; 100 } 101 102 /* perturb initial guess */ 103 PetscCall(VecGetArray(x,&xx)); 104 for (i=0; i<n; i++) { 105 v2 = xx[i]*sperturb; 106 PetscCall(VecSetValues(x,1,&i,&v2,INSERT_VALUES)); 107 } 108 PetscCall(VecRestoreArray(x,&xx)); 109 110 PetscCall(SNESSolve(snes,NULL,x)); 111 PetscCall(SNESGetIterationNumber(snes,&it)); 112 PetscCall(PetscPrintf(PETSC_COMM_SELF,"SNES iterations = %D\n\n",it)); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Free work space. All PETSc objects should be destroyed when they 116 are no longer needed. 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 119 PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r)); 120 PetscCall(VecDestroy(&F)); PetscCall(MatDestroy(&J)); 121 PetscCall(SNESDestroy(&snes)); 122 PetscCall(PetscFinalize()); 123 return 0; 124 } 125 126 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *dummy) 127 { 128 const PetscScalar *xx; 129 PetscScalar *ff,*FF,d,d2; 130 PetscInt i,n; 131 132 PetscCall(VecGetArrayRead(x,&xx)); 133 PetscCall(VecGetArray(f,&ff)); 134 PetscCall(VecGetArray((Vec)dummy,&FF)); 135 PetscCall(VecGetSize(x,&n)); 136 d = (PetscReal)(n - 1); d2 = d*d; 137 138 if (second_order) ff[0] = d*(0.5*d*(-xx[2] + 4.*xx[1] - 3.*xx[0]) - X0DOT); 139 else ff[0] = d*(d*(xx[1] - xx[0]) - X0DOT); 140 141 for (i=1; i<n-1; i++) ff[i] = d2*(xx[i-1] - 2.*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i]; 142 143 ff[n-1] = d*d*(xx[n-1] - X1); 144 PetscCall(VecRestoreArrayRead(x,&xx)); 145 PetscCall(VecRestoreArray(f,&ff)); 146 PetscCall(VecRestoreArray((Vec)dummy,&FF)); 147 return 0; 148 } 149 150 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat prejac,void *dummy) 151 { 152 const PetscScalar *xx; 153 PetscScalar A[3],d,d2; 154 PetscInt i,n,j[3]; 155 156 PetscCall(VecGetSize(x,&n)); 157 PetscCall(VecGetArrayRead(x,&xx)); 158 d = (PetscReal)(n - 1); d2 = d*d; 159 160 i = 0; 161 if (second_order) { 162 j[0] = 0; j[1] = 1; j[2] = 2; 163 A[0] = -3.*d*d*0.5; A[1] = 4.*d*d*0.5; A[2] = -1.*d*d*0.5; 164 PetscCall(MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES)); 165 } else { 166 j[0] = 0; j[1] = 1; 167 A[0] = -d*d; A[1] = d*d; 168 PetscCall(MatSetValues(prejac,1,&i,2,j,A,INSERT_VALUES)); 169 } 170 for (i=1; i<n-1; i++) { 171 j[0] = i - 1; j[1] = i; j[2] = i + 1; 172 A[0] = d2; A[1] = -2.*d2 + 2.*xx[i]; A[2] = d2; 173 PetscCall(MatSetValues(prejac,1,&i,3,j,A,INSERT_VALUES)); 174 } 175 176 i = n-1; 177 A[0] = d*d; 178 PetscCall(MatSetValues(prejac,1,&i,1,&i,&A[0],INSERT_VALUES)); 179 180 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 181 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 182 PetscCall(MatAssemblyBegin(prejac,MAT_FINAL_ASSEMBLY)); 183 PetscCall(MatAssemblyEnd(prejac,MAT_FINAL_ASSEMBLY)); 184 185 PetscCall(VecRestoreArrayRead(x,&xx)); 186 return 0; 187 } 188 189 /*TEST 190 191 test: 192 args: -n 14 -snes_monitor_short -snes_converged_reason 193 requires: !single 194 195 test: 196 suffix: 2 197 args: -n 15 -snes_monitor_short -snes_converged_reason 198 requires: !single 199 200 test: 201 suffix: 3 202 args: -n 14 -second_order -snes_monitor_short -snes_converged_reason 203 requires: !single 204 205 TEST*/ 206