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