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