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