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