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