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