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 adding the same function 97 * PETSc will be able to ignore 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 adding the same function 102 * PETSc will be able to ignore 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 162 PetscFunctionBeginUser; 163 PetscCall(VecSet(x, pfive)); 164 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 193 /* 194 Get pointers to vector data. 195 - For default PETSc vectors, VecGetArray() returns a pointer to 196 the data array. Otherwise, the routine is implementation dependent. 197 - You MUST call VecRestoreArray() when you no longer need access to 198 the array. 199 */ 200 PetscCall(VecGetArrayRead(x, &xx)); 201 PetscCall(VecGetArray(f, &ff)); 202 PetscCall(VecGetArrayRead(g, &gg)); 203 204 /* 205 Compute function 206 */ 207 PetscCall(VecGetSize(x, &n)); 208 d = (PetscReal)(n - 1); 209 d = d * d; 210 ff[0] = xx[0]; 211 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]; 212 ff[n - 1] = xx[n - 1] - 1.0; 213 214 /* 215 Restore vectors 216 */ 217 PetscCall(VecRestoreArrayRead(x, &xx)); 218 PetscCall(VecRestoreArray(f, &ff)); 219 PetscCall(VecRestoreArrayRead(g, &gg)); 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 /* ------------------------------------------------------------------- */ 223 /* 224 FormJacobian - Evaluates Jacobian matrix. 225 226 Input Parameters: 227 . snes - the SNES context 228 . x - input vector 229 . dummy - optional user-defined context (not used here) 230 231 Output Parameters: 232 . jac - Jacobian matrix 233 . B - optionally different preconditioning matrix 234 235 */ 236 237 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 238 { 239 const PetscScalar *xx; 240 PetscScalar A[3], d; 241 PetscInt i, n, j[3]; 242 243 PetscFunctionBeginUser; 244 /* 245 Get pointer to vector data 246 */ 247 PetscCall(VecGetArrayRead(x, &xx)); 248 249 /* 250 Compute Jacobian entries and insert into matrix. 251 - Note that in this case we set all elements for a particular 252 row at once. 253 */ 254 PetscCall(VecGetSize(x, &n)); 255 d = (PetscReal)(n - 1); 256 d = d * d; 257 258 /* 259 Interior grid points 260 */ 261 for (i = 1; i < n - 1; i++) { 262 j[0] = i - 1; 263 j[1] = i; 264 j[2] = i + 1; 265 A[0] = A[2] = d; 266 A[1] = -2.0 * d + 2.0 * xx[i]; 267 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 268 } 269 270 /* 271 Boundary points 272 */ 273 i = 0; 274 A[0] = 1.0; 275 276 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 277 278 i = n - 1; 279 A[0] = 1.0; 280 281 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 282 283 /* 284 Restore vector 285 */ 286 PetscCall(VecRestoreArrayRead(x, &xx)); 287 288 /* 289 Assemble matrix 290 */ 291 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 292 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 293 if (jac != B) { 294 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 295 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 296 } 297 PetscFunctionReturn(PETSC_SUCCESS); 298 } 299 300 PetscErrorCode MySNESConvergedReasonView(SNES snes, void *ctx) 301 { 302 ReasonViewCtx *monP = (ReasonViewCtx *)ctx; 303 PetscViewer viewer = monP->viewer; 304 SNESConvergedReason reason; 305 const char *strreason; 306 307 PetscFunctionBeginUser; 308 PetscCall(SNESGetConvergedReason(snes, &reason)); 309 PetscCall(SNESGetConvergedReasonString(snes, &strreason)); 310 PetscCall(PetscViewerASCIIPrintf(viewer, "Customized SNES converged reason view\n")); 311 PetscCall(PetscViewerASCIIAddTab(viewer, 1)); 312 if (reason > 0) { 313 PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", strreason)); 314 } else if (reason <= 0) { 315 PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", strreason)); 316 } 317 PetscCall(PetscViewerASCIISubtractTab(viewer, 1)); 318 PetscFunctionReturn(PETSC_SUCCESS); 319 } 320 321 PetscErrorCode MyKSPConvergedReasonView(KSP ksp, void *ctx) 322 { 323 ReasonViewCtx *monP = (ReasonViewCtx *)ctx; 324 PetscViewer viewer = monP->viewer; 325 KSPConvergedReason reason; 326 const char *reasonstr; 327 328 PetscFunctionBeginUser; 329 PetscCall(KSPGetConvergedReason(ksp, &reason)); 330 PetscCall(KSPGetConvergedReasonString(ksp, &reasonstr)); 331 PetscCall(PetscViewerASCIIAddTab(viewer, 2)); 332 PetscCall(PetscViewerASCIIPrintf(viewer, "Customized KSP converged reason view\n")); 333 PetscCall(PetscViewerASCIIAddTab(viewer, 1)); 334 if (reason > 0) { 335 PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", reasonstr)); 336 } else if (reason <= 0) { 337 PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", reasonstr)); 338 } 339 PetscCall(PetscViewerASCIISubtractTab(viewer, 3)); 340 PetscFunctionReturn(PETSC_SUCCESS); 341 } 342 343 /*TEST 344 345 test: 346 suffix: 1 347 nsize: 1 348 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 349 350 test: 351 suffix: 2 352 nsize: 1 353 args: -ksp_converged_reason_view_cancel 354 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 355 356 test: 357 suffix: 3 358 nsize: 1 359 args: -ksp_converged_reason_view_cancel -ksp_converged_reason 360 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 361 362 test: 363 suffix: 4 364 nsize: 1 365 args: -snes_converged_reason_view_cancel 366 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 367 368 test: 369 suffix: 5 370 nsize: 1 371 args: -snes_converged_reason_view_cancel -snes_converged_reason 372 filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g" 373 374 TEST*/ 375