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 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 44 PetscCheckFalse(size != 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 45 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 46 h = 1.0/(n-1); 47 comm = PETSC_COMM_WORLD; 48 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 49 Create nonlinear solver context 50 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 51 52 CHKERRQ(SNESCreate(comm,&snes)); 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 CHKERRQ(VecCreate(comm,&x)); 61 CHKERRQ(VecSetSizes(x,PETSC_DECIDE,n)); 62 CHKERRQ(VecSetFromOptions(x)); 63 CHKERRQ(VecDuplicate(x,&r)); 64 CHKERRQ(VecDuplicate(x,&F)); 65 CHKERRQ(VecDuplicate(x,&U)); 66 67 /* 68 Set function evaluation routine and vector 69 */ 70 CHKERRQ(SNESSetFunction(snes,r,FormFunction,(void*)F)); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Create matrix data structure; set Jacobian evaluation routine 74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 75 76 CHKERRQ(MatCreate(comm,&J)); 77 CHKERRQ(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n)); 78 CHKERRQ(MatSetFromOptions(J)); 79 CHKERRQ(MatSeqAIJSetPreallocation(J,3,NULL)); 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 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian,NULL)); 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Customize nonlinear solver; set runtime options 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 98 /* 99 Set an optional user-defined reasonview routine 100 */ 101 CHKERRQ(PetscViewerASCIIGetStdout(comm,&monP.viewer)); 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 CHKERRQ(SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0)); 107 } 108 CHKERRQ(SNESGetKSP(snes,&ksp)); 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 CHKERRQ(KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0)); 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 CHKERRQ(SNESSetFromOptions(snes)); 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 CHKERRQ(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 130 v = xp*xp*xp; 131 CHKERRQ(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 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 CHKERRQ(FormInitialGuess(x)); 145 CHKERRQ(SNESSolve(snes,NULL,x)); 146 CHKERRQ(SNESGetIterationNumber(snes,&its)); 147 148 /* 149 Free work space. All PETSc objects should be destroyed when they 150 are no longer needed. 151 */ 152 CHKERRQ(VecDestroy(&x)); CHKERRQ(VecDestroy(&r)); 153 CHKERRQ(VecDestroy(&U)); CHKERRQ(VecDestroy(&F)); 154 CHKERRQ(MatDestroy(&J)); CHKERRQ(SNESDestroy(&snes)); 155 /*CHKERRQ(PetscViewerDestroy(&monP.viewer));*/ 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 PetscScalar pfive = .50; 169 CHKERRQ(VecSet(x,pfive)); 170 return 0; 171 } 172 /* ------------------------------------------------------------------- */ 173 /* 174 FormFunction - Evaluates nonlinear function, F(x). 175 176 Input Parameters: 177 . snes - the SNES context 178 . x - input vector 179 . ctx - optional user-defined context, as set by SNESSetFunction() 180 181 Output Parameter: 182 . f - function vector 183 184 Note: 185 The user-defined context can contain any application-specific data 186 needed for the function evaluation (such as various parameters, work 187 vectors, and grid information). In this program the context is just 188 a vector containing the right-hand-side of the discretized PDE. 189 */ 190 191 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 192 { 193 Vec g = (Vec)ctx; 194 const PetscScalar *xx,*gg; 195 PetscScalar *ff,d; 196 PetscInt i,n; 197 198 /* 199 Get pointers to vector data. 200 - For default PETSc vectors, VecGetArray() returns a pointer to 201 the data array. Otherwise, the routine is implementation dependent. 202 - You MUST call VecRestoreArray() when you no longer need access to 203 the array. 204 */ 205 CHKERRQ(VecGetArrayRead(x,&xx)); 206 CHKERRQ(VecGetArray(f,&ff)); 207 CHKERRQ(VecGetArrayRead(g,&gg)); 208 209 /* 210 Compute function 211 */ 212 CHKERRQ(VecGetSize(x,&n)); 213 d = (PetscReal)(n - 1); d = d*d; 214 ff[0] = xx[0]; 215 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]; 216 ff[n-1] = xx[n-1] - 1.0; 217 218 /* 219 Restore vectors 220 */ 221 CHKERRQ(VecRestoreArrayRead(x,&xx)); 222 CHKERRQ(VecRestoreArray(f,&ff)); 223 CHKERRQ(VecRestoreArrayRead(g,&gg)); 224 return 0; 225 } 226 /* ------------------------------------------------------------------- */ 227 /* 228 FormJacobian - Evaluates Jacobian matrix. 229 230 Input Parameters: 231 . snes - the SNES context 232 . x - input vector 233 . dummy - optional user-defined context (not used here) 234 235 Output Parameters: 236 . jac - Jacobian matrix 237 . B - optionally different preconditioning matrix 238 239 */ 240 241 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 242 { 243 const PetscScalar *xx; 244 PetscScalar A[3],d; 245 PetscInt i,n,j[3]; 246 247 /* 248 Get pointer to vector data 249 */ 250 CHKERRQ(VecGetArrayRead(x,&xx)); 251 252 /* 253 Compute Jacobian entries and insert into matrix. 254 - Note that in this case we set all elements for a particular 255 row at once. 256 */ 257 CHKERRQ(VecGetSize(x,&n)); 258 d = (PetscReal)(n - 1); d = d*d; 259 260 /* 261 Interior grid points 262 */ 263 for (i=1; i<n-1; i++) { 264 j[0] = i - 1; j[1] = i; j[2] = i + 1; 265 A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i]; 266 CHKERRQ(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 267 } 268 269 /* 270 Boundary points 271 */ 272 i = 0; A[0] = 1.0; 273 274 CHKERRQ(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 275 276 i = n-1; A[0] = 1.0; 277 278 CHKERRQ(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 279 280 /* 281 Restore vector 282 */ 283 CHKERRQ(VecRestoreArrayRead(x,&xx)); 284 285 /* 286 Assemble matrix 287 */ 288 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 289 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 290 if (jac != B) { 291 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 292 CHKERRQ(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 293 } 294 return 0; 295 } 296 297 PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx) 298 { 299 ReasonViewCtx *monP = (ReasonViewCtx*) ctx; 300 PetscViewer viewer = monP->viewer; 301 SNESConvergedReason reason; 302 const char *strreason; 303 304 CHKERRQ(SNESGetConvergedReason(snes,&reason)); 305 CHKERRQ(SNESGetConvergedReasonString(snes,&strreason)); 306 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n")); 307 CHKERRQ(PetscViewerASCIIAddTab(viewer,1)); 308 if (reason > 0) { 309 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason)); 310 } else if (reason <= 0) { 311 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason)); 312 } 313 CHKERRQ(PetscViewerASCIISubtractTab(viewer,1)); 314 return 0; 315 } 316 317 PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx) 318 { 319 ReasonViewCtx *monP = (ReasonViewCtx*) ctx; 320 PetscViewer viewer = monP->viewer; 321 KSPConvergedReason reason; 322 const char *reasonstr; 323 324 CHKERRQ(KSPGetConvergedReason(ksp,&reason)); 325 CHKERRQ(KSPGetConvergedReasonString(ksp,&reasonstr)); 326 CHKERRQ(PetscViewerASCIIAddTab(viewer,2)); 327 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n")); 328 CHKERRQ(PetscViewerASCIIAddTab(viewer,1)); 329 if (reason > 0) { 330 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr)); 331 } else if (reason <= 0) { 332 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr)); 333 } 334 CHKERRQ(PetscViewerASCIISubtractTab(viewer,3)); 335 return 0; 336 } 337 338 /*TEST 339 340 test: 341 suffix: 1 342 nsize: 1 343 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 344 345 test: 346 suffix: 2 347 nsize: 1 348 args: -ksp_converged_reason_view_cancel 349 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 350 351 test: 352 suffix: 3 353 nsize: 1 354 args: -ksp_converged_reason_view_cancel -ksp_converged_reason 355 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 356 357 test: 358 suffix: 4 359 nsize: 1 360 args: -snes_converged_reason_view_cancel 361 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 362 363 test: 364 suffix: 5 365 nsize: 1 366 args: -snes_converged_reason_view_cancel -snes_converged_reason 367 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 368 369 TEST*/ 370