1 2 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\ 3 This example employs a user-defined reasonview routine.\n\n"; 4 5 /*T 6 Concepts: SNES^basic uniprocessor example 7 Concepts: SNES^setting a user-defined reasonview routine 8 Processors: 1 9 T*/ 10 11 #include <petscsnes.h> 12 13 /* 14 User-defined routines 15 */ 16 PETSC_EXTERN PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 17 PETSC_EXTERN PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 18 PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec); 19 PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES,void*); 20 PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP,void*); 21 22 /* 23 User-defined context for monitoring 24 */ 25 typedef struct { 26 PetscViewer viewer; 27 } ReasonViewCtx; 28 29 int main(int argc,char **argv) 30 { 31 SNES snes; /* SNES context */ 32 KSP ksp; /* KSP context */ 33 Vec x,r,F,U; /* vectors */ 34 Mat J; /* Jacobian matrix */ 35 ReasonViewCtx monP; /* monitoring context */ 36 PetscErrorCode ierr; 37 PetscInt its,n = 5,i; 38 PetscMPIInt size; 39 PetscScalar h,xp,v; 40 MPI_Comm comm; 41 42 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 43 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 44 if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 45 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 46 h = 1.0/(n-1); 47 comm = PETSC_COMM_WORLD; 48 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 49 Create nonlinear solver context 50 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 51 52 ierr = SNESCreate(comm,&snes);CHKERRQ(ierr); 53 54 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 55 Create vector data structures; set function evaluation routine 56 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 57 /* 58 Note that we form 1 vector from scratch and then duplicate as needed. 59 */ 60 ierr = VecCreate(comm,&x);CHKERRQ(ierr); 61 ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr); 62 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 63 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 64 ierr = VecDuplicate(x,&F);CHKERRQ(ierr); 65 ierr = VecDuplicate(x,&U);CHKERRQ(ierr); 66 67 /* 68 Set function evaluation routine and vector 69 */ 70 ierr = SNESSetFunction(snes,r,FormFunction,(void*)F);CHKERRQ(ierr); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create matrix data structure; set Jacobian evaluation routine 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 76 ierr = MatCreate(comm,&J);CHKERRQ(ierr); 77 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr); 78 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 79 ierr = MatSeqAIJSetPreallocation(J,3,NULL);CHKERRQ(ierr); 80 81 /* 82 Set Jacobian matrix data structure and default Jacobian evaluation 83 routine. User can override with: 84 -snes_fd : default finite differencing approximation of Jacobian 85 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 86 (unless user explicitly sets preconditioner) 87 -snes_mf_operator : form preconditioning matrix as set by the user, 88 but use matrix-free approx for Jacobian-vector 89 products within Newton-Krylov method 90 */ 91 92 ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr); 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Customize nonlinear solver; set runtime options 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 98 /* 99 Set an optional user-defined reasonview routine 100 */ 101 ierr = PetscViewerASCIIGetStdout(comm,&monP.viewer);CHKERRQ(ierr); 102 /* Just make sure we can not repeat addding the same function 103 * PETSc will be able to igore the repeated function 104 */ 105 for (i=0; i<4; i++) { 106 ierr = SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0);CHKERRQ(ierr); 107 } 108 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 109 /* Just make sure we can not repeat addding the same function 110 * PETSc will be able to igore the repeated function 111 */ 112 for (i=0; i<4; i++) { 113 ierr = KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0);CHKERRQ(ierr); 114 } 115 /* 116 Set SNES/KSP/KSP/PC runtime options, e.g., 117 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 118 */ 119 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Initialize application: 123 Store right-hand-side of PDE and exact solution 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 126 xp = 0.0; 127 for (i=0; i<n; i++) { 128 v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */ 129 ierr = VecSetValues(F,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 130 v = xp*xp*xp; 131 ierr = VecSetValues(U,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 132 xp += h; 133 } 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Evaluate initial guess; then solve nonlinear system 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 /* 139 Note: The user should initialize the vector, x, with the initial guess 140 for the nonlinear solver prior to calling SNESSolve(). In particular, 141 to employ an initial guess of zero, the user should explicitly set 142 this vector to zero by calling VecSet(). 143 */ 144 ierr = FormInitialGuess(x);CHKERRQ(ierr); 145 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 146 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 147 148 /* 149 Free work space. All PETSc objects should be destroyed when they 150 are no longer needed. 151 */ 152 ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); 153 ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = VecDestroy(&F);CHKERRQ(ierr); 154 ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); 155 /*ierr = PetscViewerDestroy(&monP.viewer);CHKERRQ(ierr);*/ 156 ierr = PetscFinalize(); 157 return ierr; 158 } 159 /* ------------------------------------------------------------------- */ 160 /* 161 FormInitialGuess - Computes initial guess. 162 163 Input/Output Parameter: 164 . x - the solution vector 165 */ 166 PetscErrorCode FormInitialGuess(Vec x) 167 { 168 PetscErrorCode ierr; 169 PetscScalar pfive = .50; 170 ierr = VecSet(x,pfive);CHKERRQ(ierr); 171 return 0; 172 } 173 /* ------------------------------------------------------------------- */ 174 /* 175 FormFunction - Evaluates nonlinear function, F(x). 176 177 Input Parameters: 178 . snes - the SNES context 179 . x - input vector 180 . ctx - optional user-defined context, as set by SNESSetFunction() 181 182 Output Parameter: 183 . f - function vector 184 185 Note: 186 The user-defined context can contain any application-specific data 187 needed for the function evaluation (such as various parameters, work 188 vectors, and grid information). In this program the context is just 189 a vector containing the right-hand-side of the discretized PDE. 190 */ 191 192 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 193 { 194 Vec g = (Vec)ctx; 195 const PetscScalar *xx,*gg; 196 PetscScalar *ff,d; 197 PetscErrorCode ierr; 198 PetscInt i,n; 199 200 /* 201 Get pointers to vector data. 202 - For default PETSc vectors, VecGetArray() returns a pointer to 203 the data array. Otherwise, the routine is implementation dependent. 204 - You MUST call VecRestoreArray() when you no longer need access to 205 the array. 206 */ 207 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 208 ierr = VecGetArray(f,&ff);CHKERRQ(ierr); 209 ierr = VecGetArrayRead(g,&gg);CHKERRQ(ierr); 210 211 /* 212 Compute function 213 */ 214 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 215 d = (PetscReal)(n - 1); d = d*d; 216 ff[0] = xx[0]; 217 for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i]; 218 ff[n-1] = xx[n-1] - 1.0; 219 220 /* 221 Restore vectors 222 */ 223 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 224 ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr); 225 ierr = VecRestoreArrayRead(g,&gg);CHKERRQ(ierr); 226 return 0; 227 } 228 /* ------------------------------------------------------------------- */ 229 /* 230 FormJacobian - Evaluates Jacobian matrix. 231 232 Input Parameters: 233 . snes - the SNES context 234 . x - input vector 235 . dummy - optional user-defined context (not used here) 236 237 Output Parameters: 238 . jac - Jacobian matrix 239 . B - optionally different preconditioning matrix 240 241 */ 242 243 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 244 { 245 const PetscScalar *xx; 246 PetscScalar A[3],d; 247 PetscErrorCode ierr; 248 PetscInt i,n,j[3]; 249 250 /* 251 Get pointer to vector data 252 */ 253 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 254 255 /* 256 Compute Jacobian entries and insert into matrix. 257 - Note that in this case we set all elements for a particular 258 row at once. 259 */ 260 ierr = VecGetSize(x,&n);CHKERRQ(ierr); 261 d = (PetscReal)(n - 1); d = d*d; 262 263 /* 264 Interior grid points 265 */ 266 for (i=1; i<n-1; i++) { 267 j[0] = i - 1; j[1] = i; j[2] = i + 1; 268 A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i]; 269 ierr = MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr); 270 } 271 272 /* 273 Boundary points 274 */ 275 i = 0; A[0] = 1.0; 276 277 ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr); 278 279 i = n-1; A[0] = 1.0; 280 281 ierr = MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr); 282 283 /* 284 Restore vector 285 */ 286 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 287 288 /* 289 Assemble matrix 290 */ 291 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 292 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 293 if (jac != B) { 294 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 295 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 296 } 297 return 0; 298 } 299 300 PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx) 301 { 302 PetscErrorCode ierr; 303 ReasonViewCtx *monP = (ReasonViewCtx*) ctx; 304 PetscViewer viewer = monP->viewer; 305 SNESConvergedReason reason; 306 const char *strreason; 307 308 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 309 ierr = SNESGetConvergedReasonString(snes,&strreason);CHKERRQ(ierr); 310 ierr = PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");CHKERRQ(ierr); 311 ierr = PetscViewerASCIIAddTab(viewer,1);CHKERRQ(ierr); 312 if (reason > 0) { 313 ierr = PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);CHKERRQ(ierr); 314 } else if (reason <= 0) { 315 ierr = PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);CHKERRQ(ierr); 316 } 317 ierr = PetscViewerASCIISubtractTab(viewer,1);CHKERRQ(ierr); 318 return 0; 319 } 320 321 PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx) 322 { 323 PetscErrorCode ierr; 324 ReasonViewCtx *monP = (ReasonViewCtx*) ctx; 325 PetscViewer viewer = monP->viewer; 326 KSPConvergedReason reason; 327 const char *reasonstr; 328 329 ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr); 330 ierr = KSPGetConvergedReasonString(ksp,&reasonstr);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIAddTab(viewer,2);CHKERRQ(ierr); 332 ierr = PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");CHKERRQ(ierr); 333 ierr = PetscViewerASCIIAddTab(viewer,1);CHKERRQ(ierr); 334 if (reason > 0) { 335 ierr = PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);CHKERRQ(ierr); 336 } else if (reason <= 0) { 337 ierr = PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);CHKERRQ(ierr); 338 } 339 ierr = PetscViewerASCIISubtractTab(viewer,3);CHKERRQ(ierr); 340 return 0; 341 } 342 343 /*TEST 344 345 test: 346 suffix: 1 347 nsize: 1 348 349 test: 350 suffix: 2 351 nsize: 1 352 args: -ksp_converged_reason_view_cancel 353 354 test: 355 suffix: 3 356 nsize: 1 357 args: -ksp_converged_reason_view_cancel -ksp_converged_reason 358 359 test: 360 suffix: 4 361 nsize: 1 362 args: -snes_converged_reason_view_cancel 363 364 test: 365 suffix: 5 366 nsize: 1 367 args: -snes_converged_reason_view_cancel -snes_converged_reason 368 369 TEST*/ 370