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 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 PetscScalar pfive = .50; 188 PetscCall(VecSet(x, pfive)); 189 return 0; 190 } 191 /* ------------------------------------------------------------------- */ 192 /* 193 FormFunction - Evaluates nonlinear function, F(x). 194 195 Input Parameters: 196 . snes - the SNES context 197 . x - input vector 198 . ctx - optional user-defined context, as set by SNESSetFunction() 199 200 Output Parameter: 201 . f - function vector 202 203 Note: 204 The user-defined context can contain any application-specific data 205 needed for the function evaluation (such as various parameters, work 206 vectors, and grid information). In this program the context is just 207 a vector containing the right-hand-side of the discretized PDE. 208 */ 209 210 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx) { 211 Vec g = (Vec)ctx; 212 const PetscScalar *xx, *gg; 213 PetscScalar *ff, d; 214 PetscInt i, n; 215 216 /* 217 Get pointers to vector data. 218 - For default PETSc vectors, VecGetArray() returns a pointer to 219 the data array. Otherwise, the routine is implementation dependent. 220 - You MUST call VecRestoreArray() when you no longer need access to 221 the array. 222 */ 223 PetscCall(VecGetArrayRead(x, &xx)); 224 PetscCall(VecGetArray(f, &ff)); 225 PetscCall(VecGetArrayRead(g, &gg)); 226 227 /* 228 Compute function 229 */ 230 PetscCall(VecGetSize(x, &n)); 231 d = (PetscReal)(n - 1); 232 d = d * d; 233 ff[0] = xx[0]; 234 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]; 235 ff[n - 1] = xx[n - 1] - 1.0; 236 237 /* 238 Restore vectors 239 */ 240 PetscCall(VecRestoreArrayRead(x, &xx)); 241 PetscCall(VecRestoreArray(f, &ff)); 242 PetscCall(VecRestoreArrayRead(g, &gg)); 243 return 0; 244 } 245 /* ------------------------------------------------------------------- */ 246 /* 247 FormJacobian - Evaluates Jacobian matrix. 248 249 Input Parameters: 250 . snes - the SNES context 251 . x - input vector 252 . dummy - optional user-defined context (not used here) 253 254 Output Parameters: 255 . jac - Jacobian matrix 256 . B - optionally different preconditioning matrix 257 258 */ 259 260 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy) { 261 const PetscScalar *xx; 262 PetscScalar A[3], d; 263 PetscInt i, n, j[3]; 264 265 /* 266 Get pointer to vector data 267 */ 268 PetscCall(VecGetArrayRead(x, &xx)); 269 270 /* 271 Compute Jacobian entries and insert into matrix. 272 - Note that in this case we set all elements for a particular 273 row at once. 274 */ 275 PetscCall(VecGetSize(x, &n)); 276 d = (PetscReal)(n - 1); 277 d = d * d; 278 279 /* 280 Interior grid points 281 */ 282 for (i = 1; i < n - 1; i++) { 283 j[0] = i - 1; 284 j[1] = i; 285 j[2] = i + 1; 286 A[0] = A[2] = d; 287 A[1] = -2.0 * d + 2.0 * xx[i]; 288 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES)); 289 } 290 291 /* 292 Boundary points 293 */ 294 i = 0; 295 A[0] = 1.0; 296 297 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 298 299 i = n - 1; 300 A[0] = 1.0; 301 302 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES)); 303 304 /* 305 Restore vector 306 */ 307 PetscCall(VecRestoreArrayRead(x, &xx)); 308 309 /* 310 Assemble matrix 311 */ 312 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 313 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 314 if (jac != B) { 315 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 316 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 317 } 318 return 0; 319 } 320 /* ------------------------------------------------------------------- */ 321 /* 322 Monitor - User-defined monitoring routine that views the 323 current iterate with an x-window plot. 324 325 Input Parameters: 326 snes - the SNES context 327 its - iteration number 328 norm - 2-norm function value (may be estimated) 329 ctx - optional user-defined context for private data for the 330 monitor routine, as set by SNESMonitorSet() 331 332 Note: 333 See the manpage for PetscViewerDrawOpen() for useful runtime options, 334 such as -nox to deactivate all x-window output. 335 */ 336 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx) { 337 MonitorCtx *monP = (MonitorCtx *)ctx; 338 Vec x; 339 340 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm)); 341 PetscCall(SNESGetSolution(snes, &x)); 342 PetscCall(VecView(x, monP->viewer)); 343 return 0; 344 } 345 346 /*TEST 347 348 test: 349 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 350 351 test: 352 suffix: 2 353 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view 354 requires: !single 355 356 test: 357 suffix: 3 358 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 359 360 test: 361 suffix: 4 362 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view 363 requires: !single 364 365 TEST*/ 366