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++) PetscCall(SNESConvergedReasonViewSet(snes, MySNESConvergedReasonView, &monP, 0)); 100 PetscCall(SNESGetKSP(snes, &ksp)); 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++) PetscCall(KSPConvergedReasonViewSet(ksp, MyKSPConvergedReasonView, &monP, 0)); 105 /* 106 Set SNES/KSP/KSP/PC runtime options, e.g., 107 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 108 */ 109 PetscCall(SNESSetFromOptions(snes)); 110 111 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 112 Initialize application: 113 Store right-hand-side of PDE and exact solution 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 116 xp = 0.0; 117 for (i = 0; i < n; i++) { 118 v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */ 119 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES)); 120 v = xp * xp * xp; 121 PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES)); 122 xp += h; 123 } 124 125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126 Evaluate initial guess; then solve nonlinear system 127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 128 /* 129 Note: The user should initialize the vector, x, with the initial guess 130 for the nonlinear solver prior to calling SNESSolve(). In particular, 131 to employ an initial guess of zero, the user should explicitly set 132 this vector to zero by calling VecSet(). 133 */ 134 PetscCall(FormInitialGuess(x)); 135 PetscCall(SNESSolve(snes, NULL, x)); 136 PetscCall(SNESGetIterationNumber(snes, &its)); 137 138 /* 139 Free work space. All PETSc objects should be destroyed when they 140 are no longer needed. 141 */ 142 PetscCall(VecDestroy(&x)); 143 PetscCall(VecDestroy(&r)); 144 PetscCall(VecDestroy(&U)); 145 PetscCall(VecDestroy(&F)); 146 PetscCall(MatDestroy(&J)); 147 PetscCall(SNESDestroy(&snes)); 148 PetscCall(PetscFinalize()); 149 return 0; 150 } 151 152 /* 153 FormInitialGuess - Computes initial guess. 154 155 Input/Output Parameter: 156 . x - the solution vector 157 */ 158 PetscErrorCode FormInitialGuess(Vec x) 159 { 160 PetscScalar pfive = .50; 161 PetscCall(VecSet(x, pfive)); 162 return 0; 163 } 164 165 /* 166 FormFunction - Evaluates nonlinear function, F(x). 167 168 Input Parameters: 169 . snes - the SNES context 170 . x - input vector 171 . ctx - optional user-defined context, as set by SNESSetFunction() 172 173 Output Parameter: 174 . f - function vector 175 176 Note: 177 The user-defined context can contain any application-specific data 178 needed for the function evaluation (such as various parameters, work 179 vectors, and grid information). In this program the context is just 180 a vector containing the right-hand-side of the discretized PDE. 181 */ 182 183 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) 184 { 185 Vec g = (Vec)ctx; 186 const PetscScalar *xx, *gg; 187 PetscScalar *ff, d; 188 PetscInt i, n; 189 190 /* 191 Get pointers to vector data. 192 - For default PETSc vectors, VecGetArray() returns a pointer to 193 the data array. Otherwise, the routine is implementation dependent. 194 - You MUST call VecRestoreArray() when you no longer need access to 195 the array. 196 */ 197 PetscCall(VecGetArrayRead(x, &xx)); 198 PetscCall(VecGetArray(f, &ff)); 199 PetscCall(VecGetArrayRead(g, &gg)); 200 201 /* 202 Compute function 203 */ 204 PetscCall(VecGetSize(x, &n)); 205 d = (PetscReal)(n - 1); 206 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); 252 d = d * d; 253 254 /* 255 Interior grid points 256 */ 257 for (i = 1; i < n - 1; i++) { 258 j[0] = i - 1; 259 j[1] = i; 260 j[2] = i + 1; 261 A[0] = A[2] = d; 262 A[1] = -2.0 * d + 2.0 * xx[i]; 263 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 264 } 265 266 /* 267 Boundary points 268 */ 269 i = 0; 270 A[0] = 1.0; 271 272 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 273 274 i = n - 1; 275 A[0] = 1.0; 276 277 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 278 279 /* 280 Restore vector 281 */ 282 PetscCall(VecRestoreArrayRead(x, &xx)); 283 284 /* 285 Assemble matrix 286 */ 287 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 288 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 289 if (jac != B) { 290 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 291 PetscCall(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 PetscCall(SNESGetConvergedReason(snes, &reason)); 304 PetscCall(SNESGetConvergedReasonString(snes, &strreason)); 305 PetscCall(PetscViewerASCIIPrintf(viewer, "Customized SNES converged reason view\n")); 306 PetscCall(PetscViewerASCIIAddTab(viewer, 1)); 307 if (reason > 0) { 308 PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", strreason)); 309 } else if (reason <= 0) { 310 PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", strreason)); 311 } 312 PetscCall(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 PetscCall(KSPGetConvergedReason(ksp, &reason)); 324 PetscCall(KSPGetConvergedReasonString(ksp, &reasonstr)); 325 PetscCall(PetscViewerASCIIAddTab(viewer, 2)); 326 PetscCall(PetscViewerASCIIPrintf(viewer, "Customized KSP converged reason view\n")); 327 PetscCall(PetscViewerASCIIAddTab(viewer, 1)); 328 if (reason > 0) { 329 PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", reasonstr)); 330 } else if (reason <= 0) { 331 PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", reasonstr)); 332 } 333 PetscCall(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