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 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 53 PetscCheckFalse(size != 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 54 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 55 h = 1.0/(n-1); 56 57 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 58 Create nonlinear solver context 59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 60 61 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 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 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 70 CHKERRQ(VecSetSizes(x,PETSC_DECIDE,n)); 71 CHKERRQ(VecSetFromOptions(x)); 72 CHKERRQ(VecDuplicate(x,&r)); 73 CHKERRQ(VecDuplicate(x,&F)); 74 CHKERRQ(VecDuplicate(x,&U)); 75 76 /* 77 Set function evaluation routine and vector 78 */ 79 CHKERRQ(SNESSetFunction(snes,r,FormFunction,(void*)F)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Create matrix data structure; set Jacobian evaluation routine 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 85 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&J)); 86 CHKERRQ(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n)); 87 CHKERRQ(MatSetFromOptions(J)); 88 CHKERRQ(MatSeqAIJSetPreallocation(J,3,NULL)); 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 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian,NULL)); 102 103 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 104 Customize nonlinear solver; set runtime options 105 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 106 107 /* 108 Set an optional user-defined monitoring routine 109 */ 110 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer)); 111 CHKERRQ(SNESMonitorSet(snes,Monitor,&monP,0)); 112 113 /* 114 Set names for some vectors to facilitate monitoring (optional) 115 */ 116 CHKERRQ(PetscObjectSetName((PetscObject)x,"Approximate Solution")); 117 CHKERRQ(PetscObjectSetName((PetscObject)U,"Exact Solution")); 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 CHKERRQ(SNESSetFromOptions(snes)); 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 CHKERRQ(SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf)); 131 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf)); 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 CHKERRQ(VecSetValues(F,1,&i,&v,INSERT_VALUES)); 142 v = xp*xp*xp; 143 CHKERRQ(VecSetValues(U,1,&i,&v,INSERT_VALUES)); 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 CHKERRQ(FormInitialGuess(x)); 157 CHKERRQ(SNESSolve(snes,NULL,x)); 158 CHKERRQ(SNESGetIterationNumber(snes,&its)); 159 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its)); 160 161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 Check solution and clean up 163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 164 165 /* 166 Check the error 167 */ 168 CHKERRQ(VecAXPY(x,none,U)); 169 CHKERRQ(VecNorm(x,NORM_2,&norm)); 170 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its)); 171 172 /* 173 Free work space. All PETSc objects should be destroyed when they 174 are no longer needed. 175 */ 176 CHKERRQ(VecDestroy(&x)); CHKERRQ(VecDestroy(&r)); 177 CHKERRQ(VecDestroy(&U)); CHKERRQ(VecDestroy(&F)); 178 CHKERRQ(MatDestroy(&J)); CHKERRQ(SNESDestroy(&snes)); 179 CHKERRQ(PetscViewerDestroy(&monP.viewer)); 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 PetscScalar pfive = .50; 193 CHKERRQ(VecSet(x,pfive)); 194 return 0; 195 } 196 /* ------------------------------------------------------------------- */ 197 /* 198 FormFunction - Evaluates nonlinear function, F(x). 199 200 Input Parameters: 201 . snes - the SNES context 202 . x - input vector 203 . ctx - optional user-defined context, as set by SNESSetFunction() 204 205 Output Parameter: 206 . f - function vector 207 208 Note: 209 The user-defined context can contain any application-specific data 210 needed for the function evaluation (such as various parameters, work 211 vectors, and grid information). In this program the context is just 212 a vector containing the right-hand-side of the discretized PDE. 213 */ 214 215 PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx) 216 { 217 Vec g = (Vec)ctx; 218 const PetscScalar *xx,*gg; 219 PetscScalar *ff,d; 220 PetscInt i,n; 221 222 /* 223 Get pointers to vector data. 224 - For default PETSc vectors, VecGetArray() returns a pointer to 225 the data array. Otherwise, the routine is implementation dependent. 226 - You MUST call VecRestoreArray() when you no longer need access to 227 the array. 228 */ 229 CHKERRQ(VecGetArrayRead(x,&xx)); 230 CHKERRQ(VecGetArray(f,&ff)); 231 CHKERRQ(VecGetArrayRead(g,&gg)); 232 233 /* 234 Compute function 235 */ 236 CHKERRQ(VecGetSize(x,&n)); 237 d = (PetscReal)(n - 1); d = d*d; 238 ff[0] = xx[0]; 239 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]; 240 ff[n-1] = xx[n-1] - 1.0; 241 242 /* 243 Restore vectors 244 */ 245 CHKERRQ(VecRestoreArrayRead(x,&xx)); 246 CHKERRQ(VecRestoreArray(f,&ff)); 247 CHKERRQ(VecRestoreArrayRead(g,&gg)); 248 return 0; 249 } 250 /* ------------------------------------------------------------------- */ 251 /* 252 FormJacobian - Evaluates Jacobian matrix. 253 254 Input Parameters: 255 . snes - the SNES context 256 . x - input vector 257 . dummy - optional user-defined context (not used here) 258 259 Output Parameters: 260 . jac - Jacobian matrix 261 . B - optionally different preconditioning matrix 262 263 */ 264 265 PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy) 266 { 267 const PetscScalar *xx; 268 PetscScalar A[3],d; 269 PetscInt i,n,j[3]; 270 271 /* 272 Get pointer to vector data 273 */ 274 CHKERRQ(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 CHKERRQ(VecGetSize(x,&n)); 282 d = (PetscReal)(n - 1); d = d*d; 283 284 /* 285 Interior grid points 286 */ 287 for (i=1; i<n-1; i++) { 288 j[0] = i - 1; j[1] = i; j[2] = i + 1; 289 A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i]; 290 CHKERRQ(MatSetValues(B,1,&i,3,j,A,INSERT_VALUES)); 291 } 292 293 /* 294 Boundary points 295 */ 296 i = 0; A[0] = 1.0; 297 298 CHKERRQ(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 299 300 i = n-1; A[0] = 1.0; 301 302 CHKERRQ(MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES)); 303 304 /* 305 Restore vector 306 */ 307 CHKERRQ(VecRestoreArrayRead(x,&xx)); 308 309 /* 310 Assemble matrix 311 */ 312 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 313 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 314 if (jac != B) { 315 CHKERRQ(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 316 CHKERRQ(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 { 338 MonitorCtx *monP = (MonitorCtx*) ctx; 339 Vec x; 340 341 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm)); 342 CHKERRQ(SNESGetSolution(snes,&x)); 343 CHKERRQ(VecView(x,monP->viewer)); 344 return 0; 345 } 346 347 /*TEST 348 349 test: 350 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always 351 352 test: 353 suffix: 2 354 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view 355 requires: !single 356 357 test: 358 suffix: 3 359 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 360 361 test: 362 suffix: 4 363 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view 364 requires: !single 365 366 TEST*/ 367