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 SNES snes; /* SNES context */ 25 KSP ksp; /* KSP context */ 26 Vec x, r, F, U; /* vectors */ 27 Mat J; /* Jacobian matrix */ 28 ReasonViewCtx monP; /* monitoring context */ 29 PetscInt its, n = 5, i; 30 PetscMPIInt size; 31 PetscScalar h, xp, v; 32 MPI_Comm comm; 33 34 PetscFunctionBeginUser; 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++) PetscCall(SNESConvergedReasonViewSet(snes, MySNESConvergedReasonView, &monP, 0)); 99 PetscCall(SNESGetKSP(snes, &ksp)); 100 /* Just make sure we can not repeat addding the same function 101 * PETSc will be able to igore the repeated function 102 */ 103 for (i = 0; i < 4; i++) PetscCall(KSPConvergedReasonViewSet(ksp, MyKSPConvergedReasonView, &monP, 0)); 104 /* 105 Set SNES/KSP/KSP/PC runtime options, e.g., 106 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 107 */ 108 PetscCall(SNESSetFromOptions(snes)); 109 110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 111 Initialize application: 112 Store right-hand-side of PDE and exact solution 113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 114 115 xp = 0.0; 116 for (i = 0; i < n; i++) { 117 v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */ 118 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES)); 119 v = xp * xp * xp; 120 PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES)); 121 xp += h; 122 } 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Evaluate initial guess; then solve nonlinear system 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 /* 128 Note: The user should initialize the vector, x, with the initial guess 129 for the nonlinear solver prior to calling SNESSolve(). In particular, 130 to employ an initial guess of zero, the user should explicitly set 131 this vector to zero by calling VecSet(). 132 */ 133 PetscCall(FormInitialGuess(x)); 134 PetscCall(SNESSolve(snes, NULL, x)); 135 PetscCall(SNESGetIterationNumber(snes, &its)); 136 137 /* 138 Free work space. All PETSc objects should be destroyed when they 139 are no longer needed. 140 */ 141 PetscCall(VecDestroy(&x)); 142 PetscCall(VecDestroy(&r)); 143 PetscCall(VecDestroy(&U)); 144 PetscCall(VecDestroy(&F)); 145 PetscCall(MatDestroy(&J)); 146 PetscCall(SNESDestroy(&snes)); 147 PetscCall(PetscFinalize()); 148 return 0; 149 } 150 151 /* 152 FormInitialGuess - Computes initial guess. 153 154 Input/Output Parameter: 155 . x - the solution vector 156 */ 157 PetscErrorCode FormInitialGuess(Vec x) { 158 PetscScalar pfive = .50; 159 PetscCall(VecSet(x, pfive)); 160 return 0; 161 } 162 163 /* 164 FormFunction - Evaluates nonlinear function, F(x). 165 166 Input Parameters: 167 . snes - the SNES context 168 . x - input vector 169 . ctx - optional user-defined context, as set by SNESSetFunction() 170 171 Output Parameter: 172 . f - function vector 173 174 Note: 175 The user-defined context can contain any application-specific data 176 needed for the function evaluation (such as various parameters, work 177 vectors, and grid information). In this program the context is just 178 a vector containing the right-hand-side of the discretized PDE. 179 */ 180 181 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) { 182 Vec g = (Vec)ctx; 183 const PetscScalar *xx, *gg; 184 PetscScalar *ff, d; 185 PetscInt i, n; 186 187 /* 188 Get pointers to vector data. 189 - For default PETSc vectors, VecGetArray() returns a pointer to 190 the data array. Otherwise, the routine is implementation dependent. 191 - You MUST call VecRestoreArray() when you no longer need access to 192 the array. 193 */ 194 PetscCall(VecGetArrayRead(x, &xx)); 195 PetscCall(VecGetArray(f, &ff)); 196 PetscCall(VecGetArrayRead(g, &gg)); 197 198 /* 199 Compute function 200 */ 201 PetscCall(VecGetSize(x, &n)); 202 d = (PetscReal)(n - 1); 203 d = d * d; 204 ff[0] = xx[0]; 205 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]; 206 ff[n - 1] = xx[n - 1] - 1.0; 207 208 /* 209 Restore vectors 210 */ 211 PetscCall(VecRestoreArrayRead(x, &xx)); 212 PetscCall(VecRestoreArray(f, &ff)); 213 PetscCall(VecRestoreArrayRead(g, &gg)); 214 return 0; 215 } 216 /* ------------------------------------------------------------------- */ 217 /* 218 FormJacobian - Evaluates Jacobian matrix. 219 220 Input Parameters: 221 . snes - the SNES context 222 . x - input vector 223 . dummy - optional user-defined context (not used here) 224 225 Output Parameters: 226 . jac - Jacobian matrix 227 . B - optionally different preconditioning matrix 228 229 */ 230 231 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) { 232 const PetscScalar *xx; 233 PetscScalar A[3], d; 234 PetscInt i, n, j[3]; 235 236 /* 237 Get pointer to vector data 238 */ 239 PetscCall(VecGetArrayRead(x, &xx)); 240 241 /* 242 Compute Jacobian entries and insert into matrix. 243 - Note that in this case we set all elements for a particular 244 row at once. 245 */ 246 PetscCall(VecGetSize(x, &n)); 247 d = (PetscReal)(n - 1); 248 d = d * d; 249 250 /* 251 Interior grid points 252 */ 253 for (i = 1; i < n - 1; i++) { 254 j[0] = i - 1; 255 j[1] = i; 256 j[2] = i + 1; 257 A[0] = A[2] = d; 258 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; 266 A[0] = 1.0; 267 268 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 269 270 i = n - 1; 271 A[0] = 1.0; 272 273 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 274 275 /* 276 Restore vector 277 */ 278 PetscCall(VecRestoreArrayRead(x, &xx)); 279 280 /* 281 Assemble matrix 282 */ 283 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 284 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 285 if (jac != B) { 286 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 287 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 288 } 289 return 0; 290 } 291 292 PetscErrorCode MySNESConvergedReasonView(SNES snes, void *ctx) { 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 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