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