1 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\ 2 This example employs a user-defined monitoring routine.\n\n"; 3 4 /* 5 Include "petscdraw.h" so that we can use PETSc drawing routines. 6 Include "petscsnes.h" so that we can use SNES solvers. Note that this 7 file automatically includes: 8 petscsys.h - base PETSc routines petscvec.h - vectors 9 petscmat.h - matrices 10 petscis.h - index sets petscksp.h - Krylov subspace methods 11 petscviewer.h - viewers petscpc.h - preconditioners 12 petscksp.h - linear solvers 13 */ 14 15 #include <petscsnes.h> 16 17 /* 18 User-defined routines 19 */ 20 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *); 21 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *); 22 extern PetscErrorCode FormInitialGuess(Vec); 23 extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *); 24 25 /* 26 User-defined context for monitoring 27 */ 28 typedef struct { 29 PetscViewer viewer; 30 } MonitorCtx; 31 32 int main(int argc, char **argv) 33 { 34 SNES snes; /* SNES context */ 35 Vec x, r, F, U; /* vectors */ 36 Mat J; /* Jacobian matrix */ 37 MonitorCtx monP; /* monitoring context */ 38 PetscInt its, n = 5, i, maxit, maxf; 39 PetscMPIInt size; 40 PetscScalar h, xp, v, none = -1.0; 41 PetscReal abstol, rtol, stol, norm; 42 43 PetscFunctionBeginUser; 44 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 45 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 46 PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 47 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 48 h = 1.0 / (n - 1); 49 50 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 51 Create nonlinear solver context 52 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 53 54 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 55 56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 57 Create vector data structures; set function evaluation routine 58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 59 /* 60 Note that we form 1 vector from scratch and then duplicate as needed. 61 */ 62 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 63 PetscCall(VecSetSizes(x, PETSC_DECIDE, n)); 64 PetscCall(VecSetFromOptions(x)); 65 PetscCall(VecDuplicate(x, &r)); 66 PetscCall(VecDuplicate(x, &F)); 67 PetscCall(VecDuplicate(x, &U)); 68 69 /* 70 Set function evaluation routine and vector 71 */ 72 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F)); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Create matrix data structure; set Jacobian evaluation routine 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 78 PetscCall(MatCreate(PETSC_COMM_WORLD, &J)); 79 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n)); 80 PetscCall(MatSetFromOptions(J)); 81 PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL)); 82 83 /* 84 Set Jacobian matrix data structure and default Jacobian evaluation 85 routine. User can override with: 86 -snes_fd : default finite differencing approximation of Jacobian 87 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 88 (unless user explicitly sets preconditioner) 89 -snes_mf_operator : form preconditioning matrix as set by the user, 90 but use matrix-free approx for Jacobian-vector 91 products within Newton-Krylov method 92 */ 93 94 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL)); 95 96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 97 Customize nonlinear solver; set runtime options 98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 99 100 /* 101 Set an optional user-defined monitoring routine 102 */ 103 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer)); 104 PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0)); 105 106 /* 107 Set names for some vectors to facilitate monitoring (optional) 108 */ 109 PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution")); 110 PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution")); 111 112 /* 113 Set SNES/KSP/KSP/PC runtime options, e.g., 114 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc> 115 */ 116 PetscCall(SNESSetFromOptions(snes)); 117 118 /* 119 Print parameters used for convergence testing (optional) ... just 120 to demonstrate this routine; this information is also printed with 121 the option -snes_view 122 */ 123 PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf)); 124 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf)); 125 126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 127 Initialize application: 128 Store right-hand-side of PDE and exact solution 129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 130 131 xp = 0.0; 132 for (i = 0; i < n; i++) { 133 v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */ 134 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES)); 135 v = xp * xp * xp; 136 PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES)); 137 xp += h; 138 } 139 140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 Evaluate initial guess; then solve nonlinear system 142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 143 /* 144 Note: The user should initialize the vector, x, with the initial guess 145 for the nonlinear solver prior to calling SNESSolve(). In particular, 146 to employ an initial guess of zero, the user should explicitly set 147 this vector to zero by calling VecSet(). 148 */ 149 PetscCall(FormInitialGuess(x)); 150 PetscCall(SNESSolve(snes, NULL, x)); 151 PetscCall(SNESGetIterationNumber(snes, &its)); 152 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its)); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Check solution and clean up 156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 157 158 /* 159 Check the error 160 */ 161 PetscCall(VecAXPY(x, none, U)); 162 PetscCall(VecNorm(x, NORM_2, &norm)); 163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its)); 164 165 /* 166 Free work space. All PETSc objects should be destroyed when they 167 are no longer needed. 168 */ 169 PetscCall(VecDestroy(&x)); 170 PetscCall(VecDestroy(&r)); 171 PetscCall(VecDestroy(&U)); 172 PetscCall(VecDestroy(&F)); 173 PetscCall(MatDestroy(&J)); 174 PetscCall(SNESDestroy(&snes)); 175 PetscCall(PetscViewerDestroy(&monP.viewer)); 176 PetscCall(PetscFinalize()); 177 return 0; 178 } 179 /* ------------------------------------------------------------------- */ 180 /* 181 FormInitialGuess - Computes initial guess. 182 183 Input/Output Parameter: 184 . x - the solution vector 185 */ 186 PetscErrorCode FormInitialGuess(Vec x) 187 { 188 PetscScalar pfive = .50; 189 PetscFunctionBeginUser; 190 PetscCall(VecSet(x, pfive)); 191 PetscFunctionReturn(PETSC_SUCCESS); 192 } 193 /* ------------------------------------------------------------------- */ 194 /* 195 FormFunction - Evaluates nonlinear function, F(x). 196 197 Input Parameters: 198 . snes - the SNES context 199 . x - input vector 200 . ctx - optional user-defined context, as set by SNESSetFunction() 201 202 Output Parameter: 203 . f - function vector 204 205 Note: 206 The user-defined context can contain any application-specific data 207 needed for the function evaluation (such as various parameters, work 208 vectors, and grid information). In this program the context is just 209 a vector containing the right-hand-side of the discretized PDE. 210 */ 211 212 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) 213 { 214 Vec g = (Vec)ctx; 215 const PetscScalar *xx, *gg; 216 PetscScalar *ff, d; 217 PetscInt i, n; 218 219 PetscFunctionBeginUser; 220 /* 221 Get pointers to vector data. 222 - For default PETSc vectors, VecGetArray() returns a pointer to 223 the data array. Otherwise, the routine is implementation dependent. 224 - You MUST call VecRestoreArray() when you no longer need access to 225 the array. 226 */ 227 PetscCall(VecGetArrayRead(x, &xx)); 228 PetscCall(VecGetArray(f, &ff)); 229 PetscCall(VecGetArrayRead(g, &gg)); 230 231 /* 232 Compute function 233 */ 234 PetscCall(VecGetSize(x, &n)); 235 d = (PetscReal)(n - 1); 236 d = d * d; 237 ff[0] = xx[0]; 238 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]; 239 ff[n - 1] = xx[n - 1] - 1.0; 240 241 /* 242 Restore vectors 243 */ 244 PetscCall(VecRestoreArrayRead(x, &xx)); 245 PetscCall(VecRestoreArray(f, &ff)); 246 PetscCall(VecRestoreArrayRead(g, &gg)); 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 /* ------------------------------------------------------------------- */ 250 /* 251 FormJacobian - Evaluates Jacobian matrix. 252 253 Input Parameters: 254 . snes - the SNES context 255 . x - input vector 256 . dummy - optional user-defined context (not used here) 257 258 Output Parameters: 259 . jac - Jacobian matrix 260 . B - optionally different preconditioning matrix 261 262 */ 263 264 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) 265 { 266 const PetscScalar *xx; 267 PetscScalar A[3], d; 268 PetscInt i, n, j[3]; 269 270 PetscFunctionBeginUser; 271 /* 272 Get pointer to vector data 273 */ 274 PetscCall(VecGetArrayRead(x, &xx)); 275 276 /* 277 Compute Jacobian entries and insert into matrix. 278 - Note that in this case we set all elements for a particular 279 row at once. 280 */ 281 PetscCall(VecGetSize(x, &n)); 282 d = (PetscReal)(n - 1); 283 d = d * d; 284 285 /* 286 Interior grid points 287 */ 288 for (i = 1; i < n - 1; i++) { 289 j[0] = i - 1; 290 j[1] = i; 291 j[2] = i + 1; 292 A[0] = A[2] = d; 293 A[1] = -2.0 * d + 2.0 * xx[i]; 294 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 295 } 296 297 /* 298 Boundary points 299 */ 300 i = 0; 301 A[0] = 1.0; 302 303 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 304 305 i = n - 1; 306 A[0] = 1.0; 307 308 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 309 310 /* 311 Restore vector 312 */ 313 PetscCall(VecRestoreArrayRead(x, &xx)); 314 315 /* 316 Assemble matrix 317 */ 318 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 319 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 320 if (jac != B) { 321 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 322 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 323 } 324 PetscFunctionReturn(PETSC_SUCCESS); 325 } 326 /* ------------------------------------------------------------------- */ 327 /* 328 Monitor - User-defined monitoring routine that views the 329 current iterate with an x-window plot. 330 331 Input Parameters: 332 snes - the SNES context 333 its - iteration number 334 norm - 2-norm function value (may be estimated) 335 ctx - optional user-defined context for private data for the 336 monitor routine, as set by SNESMonitorSet() 337 338 Note: 339 See the manpage for PetscViewerDrawOpen() for useful runtime options, 340 such as -nox to deactivate all x-window output. 341 */ 342 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx) 343 { 344 MonitorCtx *monP = (MonitorCtx *)ctx; 345 Vec x; 346 SNESConvergedReason reason; 347 348 PetscFunctionBeginUser; 349 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm)); 350 PetscCall(SNESGetConvergedReason(snes, &reason)); 351 PetscCall(SNESGetSolution(snes, &x)); 352 PetscCall(VecView(x, monP->viewer)); 353 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " converged = %s\n", SNESConvergedReasons[reason])); 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 /*TEST 358 359 test: 360 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 361 362 test: 363 suffix: 2 364 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view 365 requires: !single 366 367 test: 368 suffix: 3 369 args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 370 371 test: 372 suffix: 4 373 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view 374 requires: !single 375 376 test: 377 suffix: 5 378 filter: grep -v atol | sed -e "s/CONVERGED_ITS/DIVERGED_MAX_IT/g" | sed -e "s/CONVERGED_FNORM_RELATIVE/DIVERGED_MAX_IT/g" 379 args: -nox -snes_type {{newtonls newtontr ncg ngmres qn anderson nrichardson ms ksponly ksptransposeonly vinewtonrsls vinewtonssls fas ms}} -snes_max_it 1 380 requires: !single 381 382 TEST*/ 383