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