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