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(PetscViewerDestroy(&monP.viewer));*/ 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 PetscScalar pfive = .50; 160 PetscCall(VecSet(x, pfive)); 161 return 0; 162 } 163 /* ------------------------------------------------------------------- */ 164 /* 165 FormFunction - Evaluates nonlinear function, F(x). 166 167 Input Parameters: 168 . snes - the SNES context 169 . x - input vector 170 . ctx - optional user-defined context, as set by SNESSetFunction() 171 172 Output Parameter: 173 . f - function vector 174 175 Note: 176 The user-defined context can contain any application-specific data 177 needed for the function evaluation (such as various parameters, work 178 vectors, and grid information). In this program the context is just 179 a vector containing the right-hand-side of the discretized PDE. 180 */ 181 182 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) { 183 Vec g = (Vec)ctx; 184 const PetscScalar *xx, *gg; 185 PetscScalar *ff, d; 186 PetscInt i, n; 187 188 /* 189 Get pointers to vector data. 190 - For default PETSc vectors, VecGetArray() returns a pointer to 191 the data array. Otherwise, the routine is implementation dependent. 192 - You MUST call VecRestoreArray() when you no longer need access to 193 the array. 194 */ 195 PetscCall(VecGetArrayRead(x, &xx)); 196 PetscCall(VecGetArray(f, &ff)); 197 PetscCall(VecGetArrayRead(g, &gg)); 198 199 /* 200 Compute function 201 */ 202 PetscCall(VecGetSize(x, &n)); 203 d = (PetscReal)(n - 1); 204 d = d * d; 205 ff[0] = xx[0]; 206 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]; 207 ff[n - 1] = xx[n - 1] - 1.0; 208 209 /* 210 Restore vectors 211 */ 212 PetscCall(VecRestoreArrayRead(x, &xx)); 213 PetscCall(VecRestoreArray(f, &ff)); 214 PetscCall(VecRestoreArrayRead(g, &gg)); 215 return 0; 216 } 217 /* ------------------------------------------------------------------- */ 218 /* 219 FormJacobian - Evaluates Jacobian matrix. 220 221 Input Parameters: 222 . snes - the SNES context 223 . x - input vector 224 . dummy - optional user-defined context (not used here) 225 226 Output Parameters: 227 . jac - Jacobian matrix 228 . B - optionally different preconditioning matrix 229 230 */ 231 232 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) { 233 const PetscScalar *xx; 234 PetscScalar A[3], d; 235 PetscInt i, n, j[3]; 236 237 /* 238 Get pointer to vector data 239 */ 240 PetscCall(VecGetArrayRead(x, &xx)); 241 242 /* 243 Compute Jacobian entries and insert into matrix. 244 - Note that in this case we set all elements for a particular 245 row at once. 246 */ 247 PetscCall(VecGetSize(x, &n)); 248 d = (PetscReal)(n - 1); 249 d = d * d; 250 251 /* 252 Interior grid points 253 */ 254 for (i = 1; i < n - 1; i++) { 255 j[0] = i - 1; 256 j[1] = i; 257 j[2] = i + 1; 258 A[0] = A[2] = d; 259 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; 267 A[0] = 1.0; 268 269 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 270 271 i = n - 1; 272 A[0] = 1.0; 273 274 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 275 276 /* 277 Restore vector 278 */ 279 PetscCall(VecRestoreArrayRead(x, &xx)); 280 281 /* 282 Assemble matrix 283 */ 284 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 285 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 286 if (jac != B) { 287 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 288 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 289 } 290 return 0; 291 } 292 293 PetscErrorCode MySNESConvergedReasonView(SNES snes, void *ctx) { 294 ReasonViewCtx *monP = (ReasonViewCtx *)ctx; 295 PetscViewer viewer = monP->viewer; 296 SNESConvergedReason reason; 297 const char *strreason; 298 299 PetscCall(SNESGetConvergedReason(snes, &reason)); 300 PetscCall(SNESGetConvergedReasonString(snes, &strreason)); 301 PetscCall(PetscViewerASCIIPrintf(viewer, "Customized SNES converged reason view\n")); 302 PetscCall(PetscViewerASCIIAddTab(viewer, 1)); 303 if (reason > 0) { 304 PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", strreason)); 305 } else if (reason <= 0) { 306 PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", strreason)); 307 } 308 PetscCall(PetscViewerASCIISubtractTab(viewer, 1)); 309 return 0; 310 } 311 312 PetscErrorCode MyKSPConvergedReasonView(KSP ksp, void *ctx) { 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